18 #define ansi_escape ((char)0x1b)
19 #define ansi_bold ansi_escape<<"[1m"
20 #define ansi_normal ansi_escape<<"[0m"
21 #define ansi_red ansi_escape<<"[31m"
22 #define ansi_green ansi_escape<<"[32m"
23 #define ansi_blue ansi_escape<<"[34m"
26 #define ONE_OVER_SQRT12 0.288675
27 #define ONE_OVER_12 0.08333333333333
31 pair<double,const DCDCTrackHit *>b){
32 if (a.second->wire->ring!=b.second->wire->ring)
33 return (a.second->wire->ring>b.second->wire->ring);
34 return (a.first>b.first);
37 pair<double,const DFDCPseudo *>b){
38 if (a.second->wire->layer!=b.second->wire->layer)
39 return (a.second->wire->layer>b.second->wire->layer);
40 return (a.first>b.first);
69 gPARMS->SetDefaultParameter(
"TRKFIT:MAX_DOCA",
MAX_DOCA,
"Maximum doca for associating hit with track");
71 gPARMS->SetDefaultParameter(
"TRKFIT:HS_DEBUG_LEVEL",
HS_DEBUG_LEVEL,
"Debug verbosity level for hit selector used in track fitting (0=no debug messages)");
72 gPARMS->SetDefaultParameter(
"TRKFIT:MAKE_DEBUG_TREES",
MAKE_DEBUG_TREES,
"Create a TTree with debugging info on hit selection for the FDC and CDC");
73 gPARMS->SetDefaultParameter(
"TRKFIT:MIN_HIT_PROB_CDC",
MIN_HIT_PROB_CDC,
"Minimum probability a CDC hit may have to be associated with a track to be included in list passed to fitter");
74 gPARMS->SetDefaultParameter(
"TRKFIT:MIN_HIT_PROB_FDC",
MIN_HIT_PROB_FDC,
"Minimum probability a FDC hit may have to be associated with a track to be included in list passed to fitter");
75 gPARMS->SetDefaultParameter(
"TRKFIT:MIN_FDC_SIGMA_ANODE_CANDIDATE",
MIN_FDC_SIGMA_ANODE_CANDIDATE,
"Minimum sigma used for FDC anode hits on track candidates");
76 gPARMS->SetDefaultParameter(
"TRKFIT:MIN_FDC_SIGMA_CATHODE_CANDIDATE",
MIN_FDC_SIGMA_CATHODE_CANDIDATE,
"Minimum sigma used for FDC cathode hits on track candidates");
77 gPARMS->SetDefaultParameter(
"TRKFIT:MIN_FDC_SIGMA_ANODE_WIREBASED",
MIN_FDC_SIGMA_ANODE_WIREBASED,
"Minimum sigma used for FDC anode hits on wire-based tracks");
78 gPARMS->SetDefaultParameter(
"TRKFIT:MIN_FDC_SIGMA_CATHODE_WIREBASED",
MIN_FDC_SIGMA_CATHODE_WIREBASED,
"Minimum sigma used for FDC cathode hits on wire-based tracks");
83 loop->GetJApplication()->Lock();
85 cdchitsel= (TTree*)gROOT->FindObject(
"cdchitsel");
87 cdchitsel =
new TTree(
"cdchitsel",
"CDC Hit Selector");
88 cdchitsel->Branch(
"H", &
cdchitdbg,
"fit_type/I:p/F:theta:mass:sigma:x:y:z:s:itheta02:itheta02s:itheta02s2:dist:doca:resi:chisq:prob:sig_phi:sig_lambda:sig_pt");
91 jerr<<
" !!! WARNING !!!"<<endl;
92 jerr<<
"It appears that the cdchitsel TTree is already defined."<<endl;
93 jerr<<
"This is probably means you are running with multiple threads."<<endl;
94 jerr<<
"To avoid complication, filling of the hit selector debug"<<endl;
95 jerr<<
"trees will be disabled for this thread."<<endl;
100 fdchitsel= (TTree*)gROOT->FindObject(
"fdchitsel");
102 fdchitsel =
new TTree(
"fdchitsel",
"FDC Hit Selector");
103 fdchitsel->Branch(
"H", &
fdchitdbg,
"fit_type/I:p/F:theta:mass:sigma_anode:sigma_cathode:x:y:z:s:itheta02:itheta02s:itheta02s2:dist:doca:resi:u:u_cathodes:resic:chisq:prob:sig_phi:sig_lambda:sig_pt");
106 jerr<<
" !!! WARNING !!!"<<endl;
107 jerr<<
"It appears that the fdchitsel TTree is already defined."<<endl;
108 jerr<<
"This is probably means you are running with multiple threads."<<endl;
109 jerr<<
"To avoid complication, filling of the hit selector debug"<<endl;
110 jerr<<
"trees will be disabled for this thread."<<endl;
115 loop->GetJApplication()->Unlock();
135 const vector<DTrackFitter::Extrapolation_t> &extrapolations,
const vector<const DCDCTrackHit*> &cdchits_in, vector<const DCDCTrackHit*> &cdchits_out,
int N)
const
138 vector<pair<double,const DCDCTrackHit*> >cdchits_tmp;
141 vector<const DCDCTrackHit*> cdchits_in_sorted = cdchits_in;
149 double a=-0.003*Bz*q;
150 DVector3 mom=extrapolations[extrapolations.size()/2].momentum;
153 double a_over_p=1./p_over_a;
154 double lambda=M_PI_2-mom.Theta();
155 double tanl=tan(lambda);
156 double tanl2=tanl*tanl;
157 double sinl=
sin(lambda);
158 double cosl=cos(lambda);
159 double cosl2=cosl*cosl;
160 double sinl2=sinl*sinl;
161 double pt_over_a=cosl*p_over_a;
164 double var_lambda_res=0.;
165 double var_phi_radial=0.;
166 double var_x0=0.0,var_y0=0.0;
170 bool outermost_hit=
true;
171 vector<const DCDCTrackHit*>::const_reverse_iterator iter;
172 for(iter=cdchits_in_sorted.rbegin(); iter!=cdchits_in_sorted.rend(); iter++){
179 double z0=origin.z();
180 double d2=0.,d2_old=1.e6;
183 for (
unsigned int i=0;i<extrapolations.size();i++){
184 DVector3 trackpos=extrapolations[i].position;
185 double dz=trackpos.z()-z0;
187 d2=(trackpos-wirepos).Perp2();
190 double phi=extrapolations[i].momentum.Phi();
191 double sinphi=
sin(phi);
192 double cosphi=cos(phi);
194 double my_ux=ux*sinl-cosl*cosphi;
195 double my_uy=uy*sinl-cosl*sinphi;
196 double denom=my_ux*my_ux+my_uy*my_uy;
198 double ds=((old_trackpos.x()-origin.x()-ux*dz)*my_ux
199 +(old_trackpos.y()-origin.y()-uy*dz)*my_uy)/denom;
201 double dx=old_trackpos.x()+mom.X()/mom.Mag()*ds-wirepos.x();
202 double dy=old_trackpos.y()+mom.Y()/mom.Mag()*ds-wirepos.y();
210 double var_lambda = extrapolations[i].theta2ms_sum/3.;
212 double var_phi=var_lambda*(1.+tanl2);
217 double s_4=s_sq*s_sq;
218 double sum_s_theta_ms=extrapolations[i].s_theta_ms_sum;
219 double var_k_ms=(4./3.)*sum_s_theta_ms*sum_s_theta_ms/s_4;
220 double var_k_res=var*0.0720/double(N+4)/s_4/(cosl2*cosl2);
221 var_k=var_k_ms+var_k_res;
224 var_lambda_res=12.0*var*double(N-1)/double(N*(N+1))
228 double tan_s_2R=tan(s*cosl/(2.*pt_over_a));
229 var_phi_radial=tan_s_2R*tan_s_2R*pt_over_a*pt_over_a*var_k;
234 var_lambda+=var_lambda_res;
238 var_phi+=var_phi_radial;
241 double var_pos_ms=extrapolations[i].s_theta_ms_sum/3.;
244 double doca=
sqrt(d2);
246 double resi=dist-doca;
249 double as_over_p=s*a_over_p;
250 double sin_as_over_p=
sin(as_over_p);
251 double cos_as_over_p=cos(as_over_p);
252 double one_minus_cos_as_over_p=1-cos_as_over_p;
253 double diff1=sin_as_over_p-as_over_p*cos_as_over_p;
254 double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p;
255 double dx_dk=2.*pt_over_a*pt_over_a*(sinphi*diff2-cosphi*diff1);
256 double dx_dcosl=p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
257 double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
258 double var_x=var_x0+dx_dk*dx_dk*var_k+var_pos_ms
259 +dx_dcosl*dx_dcosl*sinl*sinl*var_lambda+dx_dphi*dx_dphi*var_phi;
261 double dy_dk=-2.*pt_over_a*pt_over_a*(sinphi*diff1+cosphi*diff2);
262 double dy_dcosl=p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
263 double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
264 double var_y=var_y0+dy_dk*dy_dk*var_k+var_pos_ms
265 +dy_dcosl*dy_dcosl*sinl*sinl*var_lambda+dy_dphi*dy_dphi*var_phi;
267 double dd_dx=dx/doca;
268 double dd_dy=dy/doca;
269 double var_d=dd_dx*dd_dx*var_x+dd_dy*dd_dy*var_y;
270 double chisq=resi*resi/(var+var_d);
273 double probability = TMath::Prob(chisq, 1);
275 pair<double,const DCDCTrackHit*>myhit;
276 myhit.first=probability;
278 cdchits_tmp.push_back(myhit);
310 jerr<<
"s="<<s<<
" t=" <<hit->
tdrift <<
" doca="<<doca<<
" dist="<<dist<<
" resi="<<resi<<
" sigma="<<
" prob="<<probability<<endl;
316 old_trackpos=trackpos;
317 s=extrapolations[i].s;
329 int old_straw=1000,old_ring=1000;
330 for (
unsigned int i=0;i<cdchits_tmp.size();i++){
331 if (cdchits_tmp[i].second->wire->ring!=old_ring ||
332 abs(cdchits_tmp[i].second->wire->straw-old_straw)==1){
333 cdchits_out.push_back(cdchits_tmp[i].second);
335 old_straw=cdchits_tmp[i].second->wire->straw;
336 old_ring=cdchits_tmp[i].second->wire->ring;
346 vector<pair<double,const DCDCTrackHit*> >cdchits_tmp;
360 vector<const DCDCTrackHit*> cdchits_in_sorted = cdchits_in;
373 double Bz0=my_step->
B.z();
374 double a=-0.003*Bz0*rt->
q;
375 double p=my_step->
mom.Mag();
377 double a_over_p=1./p_over_a;
378 double lambda=M_PI_2-my_step->
mom.Theta();
379 double cosl=cos(lambda);
380 double sinl=
sin(lambda);
381 double sinl2=sinl*sinl;
382 double cosl2=cosl*cosl;
383 double tanl=tan(lambda);
384 double tanl2=tanl*tanl;
385 double pt_over_a=cosl*p_over_a;
386 double phi=my_step->
mom.Phi();
387 double cosphi=cos(phi);
388 double sinphi=
sin(phi);
390 double var_lambda_res=0.;
391 double var_phi_radial=0.;
392 double var_x0=0.0,var_y0=0.0;
394 double one_over_beta=
sqrt(1.+mass*mass/(p*p));
395 double var_pt_factor=0.016*one_over_beta/(cosl*a);
396 double var_pt_over_pt_sq=0.;
399 int old_straw=1000,old_ring=1000;
402 bool outermost_hit=
true;
403 vector<const DCDCTrackHit*>::const_reverse_iterator iter;
404 for(iter=cdchits_in_sorted.rbegin(); iter!=cdchits_in_sorted.rend(); iter++){
419 if(!isfinite(doca))
continue;
420 if(!isfinite(s))
continue;
431 double Bratio=last_step->
B.z()/Bz0;
432 double invBratio=1./Bratio;
433 pt_over_a*=invBratio;
438 double var_lambda = last_step->
itheta02/3.;
439 double var_phi=var_lambda*(1.+tanl2);
444 double var_k_over_k_sq_res=var*p_over_a*p_over_a*0.0720/double(N+4)/(s_sq*s_sq*sinl2)/cosl2;
445 double var_k_over_k_sq_ms=var_pt_factor*var_pt_factor*last_step->
invX0/s;
447 var_pt_over_pt_sq=var_k_over_k_sq_ms+var_k_over_k_sq_res;
450 var_lambda_res=12.0*var*double(N-1)/double(N*(N+1))*cosl2*cosl2/s_sq;
453 double tan_s_2R=tan(s*cosl/(2.*pt_over_a));
454 var_phi_radial=tan_s_2R*tan_s_2R*var_pt_over_pt_sq;
459 var_lambda+=var_lambda_res;
463 var_phi+=var_phi_radial;
469 double resi = dist - doca;
475 double as_over_p=s*a_over_p;
476 double sin_as_over_p=
sin(as_over_p);
477 double cos_as_over_p=cos(as_over_p);
478 double one_minus_cos_as_over_p=1-cos_as_over_p;
479 double diff1=sin_as_over_p-as_over_p*cos_as_over_p;
480 double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p;
481 double pdx_dp=pt_over_a*(cosphi*diff1-sinphi*diff2);
482 double dx_dcosl=p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
483 double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
484 double var_x=var_x0+pdx_dp*pdx_dp*var_pt_over_pt_sq+var_pos_ms
485 +dx_dcosl*dx_dcosl*sinl*sinl*var_lambda+dx_dphi*dx_dphi*var_phi;
487 double pdy_dp=pt_over_a*(sinphi*diff1+cosphi*diff2);
488 double dy_dcosl=p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
489 double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
490 double var_y=var_y0+pdy_dp*pdy_dp*var_pt_over_pt_sq+var_pos_ms
491 +dy_dcosl*dy_dcosl*sinl*sinl*var_lambda+dy_dphi*dy_dphi*var_phi;
497 double z0=origin.z();
498 DVector3 wirepos=origin+(pos.z()-z0)/uz*dir;
499 double dd_dx=(pos.x()-wirepos.x())/doca;
500 double dd_dy=(pos.y()-wirepos.y())/doca;
501 double var_d=dd_dx*dd_dx*var_x+dd_dy*dd_dy*var_y;
503 double chisq=resi*resi/(var+var_d);
506 double probability = TMath::Prob(chisq, 1);
508 pair<double,const DCDCTrackHit*>myhit;
509 myhit.first=probability;
511 cdchits_tmp.push_back(myhit);
543 jerr<<
"s="<<s<<
" t=" <<hit->
tdrift <<
" doca="<<doca<<
" dist="<<dist<<
" resi="<<resi<<
" sigma="<<
" prob="<<probability<<endl;
554 old_straw=1000,old_ring=1000;
555 for (
unsigned int i=0;i<cdchits_tmp.size();i++){
556 if (cdchits_tmp[i].second->wire->ring!=old_ring ||
557 abs(cdchits_tmp[i].second->wire->straw-old_straw)==1){
558 cdchits_out.push_back(cdchits_tmp[i].second);
560 old_straw=cdchits_tmp[i].second->wire->straw;
561 old_ring=cdchits_tmp[i].second->wire->ring;
571 vector<pair<double,const DFDCPseudo*> >fdchits_tmp;
584 vector<const DFDCPseudo*> fdchits_in_sorted = fdchits_in;
589 const double VAR_CATHODE_STRIPS=0.000225;
590 double var_cathode = VAR_CATHODE_STRIPS;
595 double Bz0=my_step->
B.z();
596 double z0=my_step->
origin.z();
597 double a=-0.003*Bz0*rt->
q;
598 double p=my_step->
mom.Mag();
601 double a_over_p=1./p_over_a;
602 double lambda=M_PI_2-my_step->
mom.Theta();
603 double cosl=cos(lambda);
604 double cosl2=cosl*cosl;
605 double sinl=
sin(lambda);
606 double sinl2=sinl*sinl;
607 double tanl=tan(lambda);
608 double tanl2=tanl*tanl;
609 double pt_over_a=cosl*p_over_a;
610 double phi=my_step->
mom.Phi();
611 double cosphi=cos(phi);
612 double sinphi=
sin(phi);
613 double var_lambda=0.,var_phi=0.,var_lambda_res=0.;
614 double var_phi_radial=0.;
617 double var_x0=0.0,var_y0=0.0;
618 double var_pt_over_pt_sq=0.;
621 bool most_downstream_hit=
true;
622 vector<const DFDCPseudo*>::const_reverse_iterator iter;
623 for(iter=fdchits_in_sorted.rbegin(); iter!=fdchits_in_sorted.rend(); iter++){
630 if(!isfinite(doca))
continue;
631 if(!isfinite(s))
continue;
639 double dz=pos.z()-z0;
647 double u_cathodes = hit->
s;
648 double resic = u - u_cathodes;
654 double pz=last_step->
mom.z();
655 double tx=last_step->
mom.x()/pz;
656 double ty=last_step->
mom.y()/pz;
657 double tu=tx*cosa-ty*sina;
658 double alpha=atan(tu);
659 double cosalpha=cos(alpha);
663 double Bz=last_step->
B.z();
664 double Bratio=Bz/Bz0;
665 double invBratio=1./Bratio;
666 pt_over_a*=invBratio;
671 double resi = dist - doca/cosalpha;
674 double probability=0.,chisq=0.;
678 double sign=(pos.x()*cosa-pos.y()*sina-hit->
w)<0?1:-1.;
679 double ucor=0.1458*Bz*(1.-0.048*last_step->
B.Perp())*sign*doca;
684 double max_deflection=0.1458*Bz*(1.-0.048*last_step->
B.Perp())*0.5;
685 var_cathode=VAR_CATHODE_STRIPS+max_deflection*max_deflection/12.;
687 double var_tot=var_anode+var_cathode;
690 var_lambda = last_step->
itheta02/3.;
691 var_phi=var_lambda*(1.+tanl2);
693 if (most_downstream_hit){
696 double var_k_over_k_sq_res=var_tot*p_over_a*p_over_a
697 *0.0720/double(N+4)/(s_sq*s_sq)/cosl2;
700 double one_over_beta=
sqrt(1.+mass*mass/p_sq);
701 double var_pt_factor=0.016*one_over_beta/(cosl*0.003*last_step->
B.z());
702 double var_k_over_k_sq_ms=var_pt_factor*var_pt_factor*last_step->
invX0/s;
704 var_pt_over_pt_sq=var_k_over_k_sq_ms+var_k_over_k_sq_res;
707 var_lambda_res=12.0*var_tot*double(N-1)/double(N*(N+1))
711 double tan_s_2R=tan(s*cosl/(2.*pt_over_a));
712 var_phi_radial=tan_s_2R*tan_s_2R*var_pt_over_pt_sq;
714 most_downstream_hit=
false;
718 var_lambda+=var_lambda_res;
722 var_phi+=var_phi_radial;
731 double as_over_p=s*a_over_p;
732 double sin_as_over_p=
sin(as_over_p);
733 double cos_as_over_p=cos(as_over_p);
734 double one_minus_cos_as_over_p=1-cos_as_over_p;
735 double diff1=sin_as_over_p-as_over_p*cos_as_over_p;
736 double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p;
737 double pdx_dp=pt_over_a*(cosphi*diff1-sinphi*diff2);
738 double dx_ds=cosl*(cosphi*cos_as_over_p-sinphi*sin_as_over_p);
739 double ds_dcosl=dz*cosl/(sinl*sinl2);
741 =p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p)
743 double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
744 double var_z0=2.*tanl2*(var_tot)*
double(2*N-1)/double(N*(N+1));
746 double var_x=var_x0+pdx_dp*pdx_dp*var_pt_over_pt_sq+var_pos_ms
747 +dx_dcosl*dx_dcosl*sinl2*var_lambda+dx_dphi*dx_dphi*var_phi
748 +dx_ds*dx_ds*var_z0/sinl2;
750 double pdy_dp=pt_over_a*(sinphi*diff1+cosphi*diff2);
751 double dy_ds=cosl*(sinphi*cos_as_over_p+cosphi*sin_as_over_p);
753 =p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p)
755 double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
756 double var_y=var_y0+pdy_dp*pdy_dp*var_pt_over_pt_sq+var_pos_ms
757 +dy_dcosl*dy_dcosl*sinl2*var_lambda+dy_dphi*dy_dphi*var_phi
758 +dy_ds*dy_ds*var_z0/sinl2;
761 double cos2a=cosa*cosa;
762 double sin2a=sina*sina;
763 double var_d=(cos2a*var_x+sin2a*var_y)/(cosalpha*cosalpha);
764 double var_u=cos2a*var_y+sin2a*var_x;
773 chisq = resi*resi/(var_d+var_anode)+resic*resic/(var_u+var_cathode);
776 probability = TMath::Prob(chisq,2);
779 pair<double,const DFDCPseudo*>myhit;
780 myhit.first=probability;
782 fdchits_tmp.push_back(myhit);
818 jerr<<
"s="<<s<<
" doca="<<doca<<
" dist="<<dist<<
" resi="<<resi<<
" resic="<<resic<<
" chisq="<<chisq<<
" prob="<<probability<<endl;
828 int old_layer=1000,old_wire=1000;
829 for (
unsigned int i=0;i<fdchits_tmp.size();i++){
830 if (fdchits_tmp[i].second->wire->layer!=old_layer ||
831 abs(fdchits_tmp[i].second->wire->wire-old_wire)==1){
832 fdchits_out.push_back(fdchits_tmp[i].second);
834 old_wire=fdchits_tmp[i].second->wire->wire;
835 old_layer=fdchits_tmp[i].second->wire->layer;
839 const vector<DTrackFitter::Extrapolation_t> &extrapolations,
const vector<const DFDCPseudo*> &fdchits_in, vector<const DFDCPseudo*> &fdchits_out,
int N)
const
842 vector<pair<double,const DFDCPseudo*> >fdchits_tmp;
851 vector<const DFDCPseudo*> fdchits_in_sorted = fdchits_in;
857 const double VAR_CATHODE_STRIPS=0.000225;
858 double var_cathode = VAR_CATHODE_STRIPS+0.15*0.15*
ONE_OVER_12;
859 double var_tot=var_anode+var_cathode;
863 double a=-0.003*Bz*q;
864 DVector3 mom=extrapolations[extrapolations.size()/2].momentum;
867 double a_over_p=1./p_over_a;
868 double lambda=M_PI_2-mom.Theta();
869 double tanl=tan(lambda);
870 double tanl2=tanl*tanl;
871 double sinl=
sin(lambda);
872 double cosl=cos(lambda);
873 double cosl2=cosl*cosl;
874 double sinl2=sinl*sinl;
875 double pt_over_a=cosl*p_over_a;
876 double z0=extrapolations[0].position.z();
878 double var_lambda=0.,var_phi=0.,var_lambda_res=0.,var_k=0.;
879 double var_phi_radial=0.;
880 double var_x0=0.0,var_y0=0.0;
881 double var_z0=2.*tanl2*(var_tot)*
double(2*N-1)/double(N*(N+1));
884 bool most_downstream_hit=
true;
885 int last_extrapolation_index=extrapolations.size()-1;
886 vector<const DFDCPseudo*>::const_reverse_iterator iter;
887 for(iter=fdchits_in_sorted.rbegin(); iter!=fdchits_in_sorted.rend(); iter++){
889 for (
int k=last_extrapolation_index;k>=0;k--){
891 DVector3 pos=extrapolations[k].position;
892 double dz=pos.z()-z0;
893 double s=extrapolations[k].s;
896 var_lambda = extrapolations[k].theta2ms_sum/3.;
900 var_phi=var_lambda*(1.+tanl2);
903 double phi=extrapolations[k].momentum.Phi();
904 double sinphi=
sin(phi);
905 double cosphi=cos(phi);
907 if (most_downstream_hit){
910 double s_4=s_sq*s_sq;
911 double sum_s_theta_ms=extrapolations[k].s_theta_ms_sum;
912 double var_k_ms=(4./3.)*sum_s_theta_ms*sum_s_theta_ms/s_4;
913 double var_k_res=var_tot*0.0720/double(N+4)/s_4/(cosl2*cosl2);
914 var_k=var_k_ms+var_k_res;
917 var_lambda_res=12.0*var_tot*double(N-1)/double(N*(N+1))
921 double tan_s_2R=tan(s*cosl/(2.*pt_over_a));
922 var_phi_radial=tan_s_2R*tan_s_2R*pt_over_a*pt_over_a*var_k;
924 most_downstream_hit=
false;
927 var_lambda+=var_lambda_res;
931 var_phi+=var_phi_radial;
934 double var_pos_ms=extrapolations[k].s_theta_ms_sum/3.;
939 double dx=hit->
xy.X()-
x;
940 double dy=hit->
xy.Y()-
y;
941 double doca=
sqrt(dx*dx+dy*dy);
948 double u=x*sina+y*cosa;
949 double u_cathodes = hit->
s;
950 double resic = u - u_cathodes;
956 double pz=extrapolations[k].momentum.z();
957 double tx=extrapolations[k].momentum.x()/pz;
958 double ty=extrapolations[k].momentum.y()/pz;
959 double tu=tx*cosa-ty*sina;
960 double alpha=atan(tu);
961 double cosalpha=cos(alpha);
964 double resi = dist - doca/cosalpha;
967 double probability=0.,chisq=0.;
970 double as_over_p=s*a_over_p;
971 double sin_as_over_p=
sin(as_over_p);
972 double cos_as_over_p=cos(as_over_p);
973 double one_minus_cos_as_over_p=1-cos_as_over_p;
974 double diff1=sin_as_over_p-as_over_p*cos_as_over_p;
975 double diff2=one_minus_cos_as_over_p-as_over_p*sin_as_over_p;
976 double dx_ds=cosl*(cosphi*cos_as_over_p-sinphi*sin_as_over_p);
977 double ds_dcosl=dz*cosl/(sinl*sinl2);
979 =p_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p)
981 double dx_dphi=-pt_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p);
982 double dx_dk=2.*pt_over_a*pt_over_a*(sinphi*diff2-cosphi*diff1);
984 double var_x=var_x0+dx_dk*dx_dk*var_k+var_pos_ms
985 +dx_dcosl*dx_dcosl*sinl2*var_lambda+dx_dphi*dx_dphi*var_phi
986 +dx_ds*dx_ds*var_z0/sinl2;
988 double dy_dk=-2.*pt_over_a*pt_over_a*(sinphi*diff1+cosphi*diff2);
989 double dy_ds=cosl*(sinphi*cos_as_over_p+cosphi*sin_as_over_p);
991 =p_over_a*(sinphi*sin_as_over_p+cosphi*one_minus_cos_as_over_p)
993 double dy_dphi=pt_over_a*(cosphi*sin_as_over_p-sinphi*one_minus_cos_as_over_p);
994 double var_y=var_y0+dy_dk*dy_dk*var_k+var_pos_ms
995 +dy_dcosl*dy_dcosl*sinl2*var_lambda+dy_dphi*dy_dphi*var_phi
996 +dy_ds*dy_ds*var_z0/sinl2;
999 double cos2a=cosa*cosa;
1000 double sin2a=sina*sina;
1001 double var_d=(cos2a*var_x+sin2a*var_y)/(cosalpha*cosalpha);
1002 double var_u=cos2a*var_y+sin2a*var_x;
1009 chisq = resi*resi/(var_d+var_anode)+resic*resic/(var_u+var_cathode);
1012 probability = TMath::Prob(chisq,2);
1015 pair<double,const DFDCPseudo*>myhit;
1016 myhit.first=probability;
1018 fdchits_tmp.push_back(myhit);
1053 jerr<<
"s="<<s<<
" doca="<<doca<<
" dist="<<dist<<
" resi="<<resi<<
" resic="<<resic<<
" chisq="<<chisq<<
" prob="<<probability<<endl;
1067 int old_layer=1000,old_wire=1000;
1068 for (
unsigned int i=0;i<fdchits_tmp.size();i++){
1069 if (fdchits_tmp[i].second->wire->layer!=old_layer ||
1070 abs(fdchits_tmp[i].second->wire->wire-old_wire)==1){
1071 fdchits_out.push_back(fdchits_tmp[i].second);
1073 old_wire=fdchits_tmp[i].second->wire->wire;
1074 old_layer=fdchits_tmp[i].second->wire->layer;
DVector2 xy
rough x,y coordinates in lab coordinate system
DTrackHitSelectorALT2(jana::JEventLoop *loop, int32_t runnumber)
const swim_step_t * GetLastSwimStep(void) const
double MIN_FDC_SIGMA_ANODE_WIREBASED
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
DVector3 GetLastDOCAPoint(void) const
const DFDCWire * wire
DFDCWire for this wire.
void GetCDCHits(fit_type_t fit_type, const DReferenceTrajectory *rt, const vector< const DCDCTrackHit * > &cdchits_in, vector< const DCDCTrackHit * > &cdchits_out, int N=20) const
class DFDCPseudo: definition for a reconstructed point in the FDC
double GetLastDistAlongWire(void) const
The DTrackHitSelector class is a base class for algorithms that will select hits from the drift chamb...
static bool DTrackHitSelector_fdchit_cmp(pair< double, const DFDCPseudo * >a, pair< double, const DFDCPseudo * >b)
virtual ~DTrackHitSelectorALT2()
void GetFDCHits(fit_type_t fit_type, const DReferenceTrajectory *rt, const vector< const DFDCPseudo * > &fdchits_in, vector< const DFDCPseudo * > &fdchits_out, int N=20) const
const DMagneticFieldMap * bfield
double DistToRT(double x, double y, double z) const
double MIN_FDC_SIGMA_CATHODE_WIREBASED
double s
< wire position computed from cathode data, assuming the avalanche occurs at the wire ...
double GetMass(void) const
static bool DTrackHitSelector_cdchit_cmp(pair< double, const DCDCTrackHit * >a, pair< double, const DCDCTrackHit * >b)
static bool DTrackHitSelector_cdchit_in_cmp(const DCDCTrackHit *a, const DCDCTrackHit *b)
double MIN_FDC_SIGMA_CATHODE_CANDIDATE
static bool DTrackHitSelector_fdchit_in_cmp(const DFDCPseudo *a, const DFDCPseudo *b)
double MIN_FDC_SIGMA_ANODE_CANDIDATE