31 gPARMS->SetDefaultParameter(
"TRKFIND:COSMICS",
COSMICS);
34 gPARMS->SetDefaultParameter(
"TRKFIND:DEBUG_HISTS",
DEBUG_HISTS);
37 gPARMS->SetDefaultParameter(
"TRKFIND:VERBOSE",
VERBOSE);
43 gPARMS->SetDefaultParameter(
"TRKFIND:CDC_MATCH_PHI",
CDC_MATCH_PHI);
49 hCDCMatch_PairD =
new TH1F(
"CDC Pair distance",
"CDC Pair distance", 100, 0.0, 20.0);
50 hCDCMatch_Axial =
new TH1F(
"CDC Match Axial",
"CDC Match Axial", 360, -M_PI_2, M_PI_2);
52 hFDCLayerRaw =
new TH1I(
"FDC Layer Raw",
"Layer of FDC hit before track finding", 24, 0.5 , 24.5);
53 hFDCLayer =
new TH1I(
"FDC Layer",
"Layer of FDC hit after segment finding", 24, 0.5 , 24.5);
54 hFDCLayerFirst =
new TH1I(
"First FDC Layer",
"First Layer of FDC segment finding", 24, 0.5 , 24.5);
76 for (
unsigned int i=0;i<
cdc_tracks.size();i++){
83 for (
unsigned int i=0;i<4;i++){
113 vector<pair<unsigned int,unsigned int> > pairs;
114 for (
unsigned int i=0;i<
axial_hits.size()-1;i++){
115 for (
unsigned int j=i+1;j<
axial_hits.size();j++){
118 int r1=first_wire->
ring;
119 int r2=second_wire->
ring;
120 int s1=first_wire->
straw;
121 int s2=second_wire->
straw;
122 double d=(first_wire->
origin-second_wire->
origin).Perp();
125 || (r1 == r2 && abs(s1-s2)==1)){
126 pair <unsigned int,unsigned int> mypair(i,j);
127 if(
VERBOSE) jout<<
" Pair formed R" << r1 <<
" S" << s1 <<
" ,R" << r2 <<
" S" << s2 << endl;
128 pairs.push_back(mypair);
134 for (
unsigned int i=0;i<pairs.size();i++){
137 vector<const DCDCTrackHit *>neighbors;
138 set<unsigned int> used;
141 neighbors.push_back(
axial_hits[pairs[i].first].hit);
142 neighbors.push_back(
axial_hits[pairs[i].second].hit);
143 used.insert(pairs[i].first);
144 used.insert(pairs[i].second);
145 for (
unsigned int j=i+1;j<pairs.size();j++){
147 unsigned int first=pairs[j].first;
148 unsigned int second=pairs[j].second;
149 set<unsigned int>::iterator it1, it2;
150 it1 = used.find(first);
151 it2 = used.find(second);
153 if (it1 != used.end() && it2 == used.end()){
158 else if (it2 != used.end() && it1 == used.end()){
175 unsigned int iHit = 1;
176 origin=neighbors[0]->wire->origin;
177 for (iHit = 1; iHit < neighbors.size(); iHit++){
178 dir+=neighbors[iHit]->wire->origin-origin;
180 if(dir.Mag() != 0.) dir.SetMag(1.);
183 for (
unsigned int iHit = 0; iHit < neighbors.size(); iHit++) dir += neighbors[iHit]->wire->origin;
184 if(dir.Mag() != 0.) dir.SetMag(1.);
188 jout <<
" Axial Segment Formed: " << endl;
189 for(
unsigned int jj = 0; jj<neighbors.size(); jj++){
190 jout <<
" R" << neighbors[jj]->wire->ring <<
" S" << neighbors[jj]->wire->straw << endl;
192 jout <<
"Dir" << endl;
209 if (num_axial<1)
return false;
210 for (
unsigned int i=0;i<num_axial;i++){
217 for (
unsigned int j=i+1;j<num_axial;j++){
220 while (dphi>M_PI) dphi-=2*M_PI;
221 while (dphi<-M_PI) dphi+=2*M_PI;
227 if ( fabs(dphi) < matchphi){
241 jout <<
" Axial track Formed: pos vhat" << endl;
242 pos0.Print(); vhat.Print();
243 for(
unsigned int jj = 0; jj<mytrack.
axial_hits.size(); jj++){
244 jout <<
" R" << mytrack.
axial_hits[jj]->wire->ring <<
" S" << mytrack.
axial_hits[jj]->wire->straw << endl;
249 bool got_match=
false;
250 for (
unsigned int j=0;j<
axial_hits.size();j++){
271 vhat.SetXYZ(0., 0., 0.);
272 for (
unsigned int iHit = 1; iHit < mytrack.
axial_hits.size(); iHit++){
294 if (
VERBOSE) jout <<
" num_axial num_stereo " << num_axial <<
" " << num_stereo << endl;
295 if (num_stereo>0 && num_stereo+num_axial>4){
313 double vhat_dot_uhat=vhat.Dot(uhat);
314 double scale=1./(1.-vhat_dot_uhat*vhat_dot_uhat);
315 double s=scale*(vhat_dot_uhat*diff.Dot(vhat)-diff.Dot(uhat));
316 double t=scale*(diff.Dot(vhat)-vhat_dot_uhat*diff.Dot(uhat));
317 if (t<0)
return false;
318 double d=(diff+s*uhat-t*vhat).Mag();
320 if (d<cut)
return true;
331 double diffCrosstdir=diff.X()*tdir.Y()-diff.Y()*tdir.X();
332 double diffCrosswdir=diff.X()*wdir.Y()-diff.Y()*wdir.X();
333 double wdirCrosstdir=wdir.X()*tdir.Y()-wdir.Y()*tdir.X();
334 double w=diffCrosstdir/wdirCrosstdir;
335 double t=diffCrosswdir/wdirCrosstdir;
336 if(t<0.)
return false;
337 if(fabs(w)<hit->
wire->
L/2.)
return true;
347 double vx=this->
dir.x();
348 double vy=this->
dir.y();
353 double sumv=0,sumx=0,sumy=0,sumz=0,sumxx=0,sumyy=0,sumxz=0,sumyz=0;
354 for (
unsigned int i=0;i<this->
stereo_hits.size();i++){
358 double ux_s=dir_s.x();
359 double uy_s=dir_s.y();
360 double dx=xa-origin_s.x();
361 double dy=ya-origin_s.y();
362 double s=(dx*vy-dy*vx)/(ux_s*vy-uy_s*vx);
364 double x=pos1.x(),
y=pos1.y(),
z=pos1.z();
366 if (
z>17.0 &&
z<167.0)
384 const double EPS=1
e-4;
385 double xdenom=sumv*sumxz-sumx*sumz;
386 if (fabs(xdenom)<EPS)
return VALUE_OUT_OF_RANGE;
388 double ydenom=sumv*sumyz-sumy*sumz;
389 if (fabs(ydenom)<EPS)
return VALUE_OUT_OF_RANGE;
391 double xtemp=sumv*sumxx-sumx*sumx;
392 double xslope=xtemp/xdenom;
393 double ytemp=sumv*sumyy-sumy*sumy;
394 double yslope=ytemp/ydenom;
397 double z0y=(sumyy*sumz-sumy*sumyz)/ytemp;
400 double delta_z=(yslope>0)?0.5:-0.5;
403 this->
z=z0y+ya/yslope+delta_z;
424 double vhat_dot_diff=diff.Dot(vhat);
425 double uhat_dot_diff=diff.Dot(uhat);
426 double uhat_dot_vhat=uhat.Dot(vhat);
427 double D=1.-uhat_dot_vhat*uhat_dot_vhat;
428 double N=uhat_dot_vhat*vhat_dot_diff-uhat_dot_diff;
429 double N1=vhat_dot_diff-uhat_dot_vhat*uhat_dot_diff;
434 if (poca!=NULL) *poca=pos1+s*uhat;
452 double vhat_dot_diff=diff.Dot(vhat);
453 double uhat_dot_diff=diff.Dot(uhat);
454 double uhat_dot_vhat=uhat.Dot(vhat);
455 double D=1.-uhat_dot_vhat*uhat_dot_vhat;
456 double N=uhat_dot_vhat*vhat_dot_diff-uhat_dot_diff;
457 double N1=vhat_dot_diff-uhat_dot_vhat*uhat_dot_diff;
462 if (poca!=NULL) *poca=pos+s*uhat;
471 if (
fdc_hits.size()==0)
return false;
472 unsigned int num_hits=
fdc_hits.size();
481 int old_layer=
fdc_hits[0].hit->wire->layer;
482 vector<unsigned int>x_list;
484 for (
unsigned int i=0;i<num_hits;i++){
488 if (
fdc_hits[i].hit->wire->layer!=old_layer){
491 old_layer=
fdc_hits[i].hit->wire->layer;
493 x_list.push_back(num_hits);
495 unsigned int start=0;
497 while (start<x_list.size()-1){
499 for (
unsigned int i=x_list[start];i<x_list[start+1];i++){
508 double z=
fdc_hits[i].hit->wire->origin.z();
511 vector<const DFDCPseudo*>neighbors;
512 neighbors.push_back(
fdc_hits[i].hit);
513 unsigned int match=0;
514 double delta,delta_min=1000.;
515 for (
unsigned int k=0;k<x_list.size()-1;k++){
519 for (
unsigned int m=x_list[k];m<x_list[k+1];m++){
521 delta=(XY-
fdc_hits[m].hit->xy).Mod();
522 double delta_z=fabs(z-
fdc_hits[m].hit->wire->origin.z());
523 if (delta<delta_min){
525 if (delta<MATCH_RADIUS && delta_z<11.0) {
535 neighbors.push_back(
fdc_hits[match].hit);
538 unsigned int num_neighbors=neighbors.size();
541 for (
unsigned int k=0;k<num_hits;k++){
543 for (
unsigned int j=0;j<num_neighbors;j++){
544 delta=(
fdc_hits[k].hit->xy-neighbors[j]->xy).Mod();
546 if (delta<ADJACENT_MATCH_RADIUS &&
547 abs(neighbors[j]->wire->wire-
fdc_hits[k].hit->wire->wire)<=1
548 && neighbors[j]->wire->origin.z()==
fdc_hits[k].hit->wire->origin.z()){
550 neighbors.push_back(
fdc_hits[k].hit);
556 if (neighbors.size()>3){
557 unsigned int packNum=(neighbors[0]->wire->layer-1)/6;
559 for (
size_t iN = 0; iN < neighbors.size(); iN++){
560 hFDCLayer->Fill(neighbors[iN]->wire->layer);
569 while (start<x_list.size()-1){
570 if (
fdc_hits[x_list[start]].used==
false)
break;
583 const double LINK_MATCH_RADIUS=7.0;
586 vector<const DFDCPseudo *>myhits;
589 for (
unsigned int i=0;i<4;i++){
592 unsigned i_plus_1=i+1;
599 for (
unsigned int k=0;k<
fdc_segments[i_plus_1].size();k++){
601 double z=
fdc_segments[i_plus_1][k].hits[0]->wire->origin.z();
604 if ((proj-
fdc_segments[i_plus_1][k].hits[0]->xy).Mod()<LINK_MATCH_RADIUS){
606 myhits.insert(myhits.end(),
fdc_segments[i_plus_1][k].hits.begin(),
609 unsigned int i_plus_2=i_plus_1+1;
616 for (
unsigned int m=0;m<
fdc_segments[i_plus_2].size();m++){
619 proj.Set(x0+tx*z,y0+ty*z);
621 if ((proj-
fdc_segments[i_plus_2][m].hits[0]->xy).Mod()<LINK_MATCH_RADIUS){
623 myhits.insert(myhits.end(),
fdc_segments[i_plus_2][m].hits.begin(),
626 unsigned int i_plus_3=i_plus_2+1;
633 for (
unsigned int n=0;n<
fdc_segments[i_plus_3].size();n++){
636 proj.Set(x0+tx*z,y0+ty*z);
638 if ((proj-
fdc_segments[i_plus_3][n].hits[0]->xy).Mod()<LINK_MATCH_RADIUS){
640 myhits.insert(myhits.end(),
fdc_segments[i_plus_3][n].hits.begin(),
658 if (myhits.size()>0){
669 for (
unsigned int i=0;i<
fdc_tracks.size()-1;i++){
671 size_t last_index_1=
fdc_tracks[i].hits.size()-1;
672 int first_pack_1=(
fdc_tracks[i].hits[0]->wire->layer-1)/6;
673 int last_pack_1=(
fdc_tracks[i].hits[last_index_1]->wire->layer-1)/6;
674 for (
unsigned int j=i+1;j<
fdc_tracks.size();j++){
675 size_t last_index_2=
fdc_tracks[j].hits.size()-1;
676 int first_pack_2=(
fdc_tracks[j].hits[0]->wire->layer-1)/6;
677 int last_pack_2=(
fdc_tracks[j].hits[last_index_2]->wire->layer-1)/6;
679 if (last_pack_1<first_pack_2 || first_pack_1 > last_pack_2){
680 double z=
fdc_tracks[j].hits[0]->wire->origin.z();
682 double diff=(
fdc_tracks[j].hits[0]->xy-proj).Mod();
684 if (diff<MATCH_RADIUS){
687 if (last_pack_1<first_pack_2){
709 for (
unsigned int i=0;i<
fdc_tracks.size();i++){
711 size_t last_index_1=
fdc_tracks[i].hits.size()-1;
712 int first_pack_1=(
fdc_tracks[i].hits[0]->wire->layer-1)/6;
713 int last_pack_1=(
fdc_tracks[i].hits[last_index_1]->wire->layer-1)/6;
714 for (
unsigned int j=0;j<4;j++){
717 int pack_2=(
fdc_segments[j][k].hits[0]->wire->layer-1)/6;
718 if (pack_2<first_pack_1 || pack_2>last_pack_1){
721 double diff=(
fdc_segments[j][k].hits[0]->xy-proj).Mod();
723 if (diff<MATCH_RADIUS){
727 if (pack_2<first_pack_1){
769 for (
unsigned int i=0;i<hits.size();i++){
770 double cosa=hits[i]->wire->udir.y();
771 double sina=hits[i]->wire->udir.x();
772 double x=hits[i]->xy.X();
773 double y=hits[i]->xy.Y();
774 double z=hits[i]->wire->origin.z();
775 double sig2x=cosa*cosa/12+sina*sina*sig2v;
776 double sig2y=sina*sina/12+cosa*cosa*sig2v;
777 double one_over_var1=1/sig2y;
778 double one_over_var2=1/sig2x;
781 S1z+=z*one_over_var1;
782 S1y+=y*one_over_var1;
783 S1zz+=z*z*one_over_var1;
784 S1zy+=z*y*one_over_var1;
787 S2z+=z*one_over_var2;
788 S2x+=x*one_over_var2;
789 S2zz+=z*z*one_over_var2;
790 S2zx+=z*x*one_over_var2;
792 double D1=S1*S1zz-S1z*S1z;
793 double y_intercept=(S1zz*S1y-S1z*S1zy)/D1;
794 double y_slope=(S1*S1zy-S1z*S1y)/D1;
795 double D2=S2*S2zz-S2z*S2z;
796 double x_intercept=(S2zz*S2x-S2z*S2zx)/D2;
797 double x_slope=(S2*S2zx-S2z*S2x)/D2;
799 return DMatrix4x1(x_intercept,y_intercept,x_slope,y_slope);
810 double dot=mydir.Dot(norm);
811 if (fabs(dot)<1
e-16)
return false;
812 double s=(origin-pos).Dot(norm)/dot;
825 double denom=dir.Mag();
826 double ux=dir.x()/denom;
827 double uy=dir.y()/denom;
828 double uz=dir.z()/denom;
831 double ux2_plus_uy2=ux2+uy2;
835 double A=ux2_plus_uy2*R*R-uy2*x0*x0-ux2*y0*y0+2.*ux*uy*x0*y0;
836 if (A<0)
return false;
838 double t0=-(x0*ux+y0*uy)/ux2_plus_uy2;
839 double dt=
sqrt(A)/ux2_plus_uy2;
841 out1.SetXYZ(x0+ux*tplus,y0+uy*tplus,z0+uz*tplus);
843 out2.SetXYZ(x0+ux*tminus,y0+uy*tminus,z0+uz*tminus);
#define ADJACENT_MATCH_RADIUS
void AddHit(const DCDCTrackHit *hit)
bool FindIntersectionWithPlane(const DVector3 &origin, const DVector3 &norm, const DVector3 &pos, const DVector3 &dir, DVector3 &outpos) const
vector< cdc_track_t > cdc_tracks
bool FindFDCSegments(void)
vector< const DCDCTrackHit * > axial_hits
const DFDCWire * wire
DFDCWire for this wire.
bool MatchCDCStereoHit(const DVector3 &tdir, const DVector3 &t0, const DCDCTrackHit *hit)
vector< cdc_hit_t > axial_hits
class DFDCPseudo: definition for a reconstructed point in the FDC
bool DTrackFinder_fdc_hit_cmp(const DTrackFinder::fdc_hit_t &a, const DTrackFinder::fdc_hit_t &b)
vector< cdc_hit_t > stereo_hits
bool LinkFDCSegments(void)
vector< fdc_segment_t > fdc_tracks
vector< fdc_segment_t > fdc_segments[4]
vector< fdc_hit_t > fdc_hits
bool DTrackFinder_cdc_hit_cosmics_cmp(const DCDCTrackHit *a, const DCDCTrackHit *b)
bool LinkCDCSegments(void)
bool MatchCDCHit(const DVector3 &vhat, const DVector3 &pos0, const DCDCTrackHit *hit, double cut)
DMatrix4x1 FindStateVector(void) const
jerror_t FindStateVector(bool IsCosmics=false)
bool FindAxialSegments(void)
double CDC_COSMIC_MATCH_PHI
bool FindIntersectionsWithCylinder(double R, const DVector3 &dir, const DVector3 &pos, DVector3 &out1, DVector3 &out2) const
vector< cdc_segment_t > axial_segments
double FindDoca(double z, const DMatrix4x1 &S, const DVector3 &wdir, const DVector3 &origin, DVector3 *poca=NULL) const
vector< const DCDCTrackHit * > stereo_hits
bool DTrackFinder_cdc_hit_cmp(const DCDCTrackHit *a, const DCDCTrackHit *b)