19 #include "JANA/JGeometry.h"
39 gPARMS->SetDefaultParameter(
"TRKFIND:MAX_SEED_DIST", MAX_SEED_DIST);
48 TARGET_Z_MIN = 65.0 - 15.0;
49 TARGET_Z_MAX = 65.0 + 15.0;
52 MAX_HIT_DIST2 = MAX_HIT_DIST*MAX_HIT_DIST;
83 vector<DFDCSeed> seeds;
87 for(
unsigned int i=0; i<seeds.size(); i++){
89 if(debug_level>3)
_DBG_<<
"----- Seed "<<i<<
" ------"<<endl;
90 if(!seed.
valid)
continue;
98 double phi = seed.
phi;
100 double theta = seed.
theta;
103 if(debug_level>3)
_DBG_<<
"p_trans="<<p_trans<<
" phi="<<phi<<
" theta="<<theta<<
" p="<<p_trans/
sin(theta)<<endl;
108 pos.SetXYZ(0.0, 0.0, z_vertex);
109 mom.SetMagThetaPhi(p_trans/
sin(theta), theta, phi);
114 _data.push_back(can);
126 for(
unsigned int i=0; i<fdctrkhits.size(); i++){
127 delete fdctrkhits[i];
132 vector<const DFDCIntersection*> fdcintersects;
133 loop->Get(fdcintersects);
136 for(
unsigned int i=0; i<fdcintersects.size(); i++){
139 trkhit->
hit = fdcintersects[i];
140 trkhit->
flags = NOISE;
142 fdctrkhits.push_back(trkhit);
147 for(
unsigned int i=0; i<fdctrkhits.size(); i++){
149 for(
unsigned int j=i+1; j<fdctrkhits.size(); j++){
151 if(!(trkhit2->
flags & NOISE))
continue;
152 double d2 = trkhit1->
Dist2(trkhit2);
153 double deltaz = fabs(trkhit1->
hit->
pos.Z() - trkhit2->
hit->
pos.Z());
154 if(debug_level>4)
_DBG_<<
" -- Dist hits "<<i<<
" and "<<j<<
" : deltaR="<<
sqrt(d2)<<
" deltaZ="<<deltaz<<endl;
155 if(d2<MAX_HIT_DIST2 && deltaz<12.0){
156 trkhit1->
flags &= ~NOISE;
157 trkhit2->
flags &= ~NOISE;
176 unsigned int last_available_hits = 0;
177 for(
int i=0; i<12; i++){
180 unsigned int Navailable_hits = NumAvailableHits();
181 if(debug_level>3)
_DBG_<<
"Navailable_hits="<<Navailable_hits<<
" last_available_hits="<<last_available_hits<<endl;
182 if(Navailable_hits<3)
break;
183 if(Navailable_hits==last_available_hits)
break;
184 last_available_hits = Navailable_hits;
187 hough.SetLimits(-400.0, +400.0, -400.0, +400.0, 100, 100);
191 for(
unsigned int i=0; i<fdctrkhits.size(); i++){
193 if(fdctrkhit->
flags&USED)
continue;
194 if(fdctrkhit->
flags&NOISE)
continue;
195 if(fdctrkhit->
flags&OUT_OF_TIME)
continue;
198 hough.AddPoint(hit->
pos.X(), hit->
pos.Y());
202 if(debug_level>3)
_DBG_<<
"Rough: GetMaxBinContent="<<hough.GetMaxBinContent()<<
" x0="<<Ro.X()<<
" y0="<<Ro.Y()<<endl;
203 if(debug_level>6)hough.PrintHist();
204 if(hough.GetMaxBinContent()<10.0)
break;
208 hough.SetLimits(Ro.X()-width, Ro.X()+width, Ro.Y()-width, Ro.Y()+width, 100, 100);
210 if(debug_level>3)
_DBG_<<
"Medium: GetMaxBinContent="<<hough.GetMaxBinContent()<<
" x0="<<Ro.X()<<
" y0="<<Ro.Y()<<endl;
211 if(debug_level>5)hough.PrintHist();
215 hough.SetLimits(Ro.X()-width, Ro.X()+width, Ro.Y()-width, Ro.Y()+width, 100, 100);
217 if(debug_level>3)
_DBG_<<
"Fine: GetMaxBinContent="<<hough.GetMaxBinContent()<<
" x0="<<Ro.X()<<
" y0="<<Ro.Y()<<endl;
218 if(debug_level>5)hough.PrintHist();
225 seed.
theta = M_PI/2.0;
236 if(seed.
hits.size()<3)
continue;
239 seed.
phi = Ro.Phi() - seed.
q*M_PI/2.0;
245 seeds.push_back(seed);
258 for(
unsigned int i=0; i<fdctrkhits.size(); i++){
260 fdctrkhit->
flags &= ~ON_CIRCLE;
261 if(fdctrkhit->
flags&USED)
continue;
262 if(fdctrkhit->
flags&NOISE)
continue;
263 if(fdctrkhit->
flags&OUT_OF_TIME)
continue;
274 double dist = fabs(g.X()*Ro_minus_h.Y() - g.Y()*Ro_minus_h.X());
278 if(debug_level>3)
_DBG_<<
"dist="<<dist<<endl;
280 fdctrkhit->
flags &= ~VALID_HIT;
284 fdctrkhit->
flags |= ON_CIRCLE | VALID_HIT;
285 seed.
hits.push_back(fdctrkhit);
295 DVector2 last_dir = -1.0*Ro/Ro.Mod();
296 double last_phi = 0.0;
300 for(
unsigned int i=0; i<seed.
hits.size(); i++){
307 p_minus_Ro/=p_minus_Ro.Mod();
308 double delta_phi = p_minus_Ro.Phi() - last_dir.Phi();
309 while(delta_phi>+M_PI)delta_phi-=2.0*M_PI;
310 while(delta_phi<-M_PI)delta_phi+=2.0*M_PI;
311 fdctrkhit->
phi_hit = last_phi + delta_phi;
313 last_dir = p_minus_Ro;
314 if(debug_level>3)
_DBG_<<
"phi_hit="<<fdctrkhit->
phi_hit<<
" z="<<fdctrkhit->
hit->
pos.Z()<<
" phi_hit/z="<<fdctrkhit->
phi_hit/(fdctrkhit->
hit->
pos.Z()-65.0)<<
" delta_phi="<<delta_phi<<endl;
316 Nq += fdctrkhit->
phi_hit>0.0 ? +1:-1;
322 for(
unsigned int i=0; i<seed.
hits.size(); i++){
325 Nq += fdctrkhit->
phi_hit>0.0 ? +1:-1;
330 if(Nq<0)seed.
q = -seed.
q;
341 for(
unsigned int i=0; i<fdctrkhits.size(); i++){
343 if(fdctrkhit->
flags&USED)
continue;
344 if(fdctrkhit->
flags&NOISE)
continue;
345 if(fdctrkhit->
flags&OUT_OF_TIME)
continue;
346 if(fdctrkhit->
flags&CANT_BE_IN_SEED)
continue;
358 FindTheta(seed, TARGET_Z_MIN, TARGET_Z_MAX);
364 if(debug_level>3)
_DBG_<<
"Seed z-vertex outside of target range (z="<<seed.
z_vertex<<
" TARGET_Z_MIN="<<TARGET_Z_MIN<<
" TARGET_Z_MAX="<<TARGET_Z_MAX<<endl;
371 for(
unsigned int i=0; i<seed.
hits.size(); i++){
373 if(trkhit->
flags&IN_Z_RANGE)trkhit->
flags |= USED;
379 for(
unsigned int i=0; i<seed.
hits.size(); i++){
381 if(trkhit->
flags&ON_CIRCLE)trkhit->
flags |= CANT_BE_IN_SEED;
404 unsigned int Nbins = 1000;
406 for(
unsigned int i=0; i<
Nbins; i++)hist[i] = 0;
407 double bin_width = 2.0*M_PI/(double)Nbins;
408 double hist_low_limit = -M_PI;
411 double &r0 = seed.
r0;
412 for(
unsigned int i=0; i<seed.
hits.size(); i++){
414 if(!trkhit->
flags&VALID_HIT)
continue;
415 if(!trkhit->
flags&ON_CIRCLE)
continue;
419 if(seed.
q<0.0)alpha = -
alpha;
420 double z_hit = trkhit->
hit->
pos.Z();
421 double tmin = atan2(alpha, z_hit - target_z_min);
422 double tmax = atan2(alpha, z_hit - target_z_max);
428 if(debug_level>3)
_DBG_<<
" -- phi_hit="<<trkhit->
phi_hit<<
" z_hit="<<z_hit<<endl;
429 if(debug_level>3)
_DBG_<<
" -- tmin="<<tmin<<
" tmax="<<tmax<<endl;
436 unsigned int imin = (
unsigned int)floor((tmin-hist_low_limit)/bin_width);
437 unsigned int imax = (
unsigned int)floor((tmax-hist_low_limit)/bin_width);
441 if(imin>=Nbins)
continue;
444 if(imin>=Nbins)imin=Nbins-1;
445 if(imax>=Nbins)imax=Nbins-1;
448 for(
unsigned int j=imin; j<=imax; j++)hist[j]++;
452 unsigned int istart=0;
454 for(
unsigned int i=1; i<
Nbins; i++){
455 if(hist[i]>hist[istart]){
457 if(debug_level>3)
_DBG_<<
" -- istart="<<istart<<
" (theta="<<hist_low_limit + bin_width*(0.5+(double)istart)<<
" , N="<<hist[i]<<
")"<<endl;
459 if(hist[i] == hist[istart])iend = i;
463 if(hist[istart]==0.0){
464 if(debug_level>3)
_DBG_<<
" - No entries in theta hist. Seed aborted."<<endl;
469 seed.
theta_min = hist_low_limit + bin_width*(0.5+(double)istart);
470 seed.
theta_max = hist_low_limit + bin_width*(0.5+(double)iend);
472 if(debug_level>3)
_DBG_<<
"istart="<<istart<<
" iend="<<iend<<
" theta_min="<<seed.
theta_min<<
" theta_max="<<seed.
theta_max<<endl;
476 for(
unsigned int i=0; i<seed.
hits.size(); i++){
478 trkhit->
flags &= ~IN_THETA_RANGE;
479 if(!trkhit->
flags&VALID_HIT)
continue;
480 if(!trkhit->
flags&ON_CIRCLE)
continue;
483 trkhit->
flags |= IN_THETA_RANGE;
494 if(seed.
phi > 2.0*M_PI)seed.
phi -= 2.0*M_PI;
515 unsigned int Nbins = 300;
517 for(
unsigned int i=0; i<
Nbins; i++)hist[i] = 0;
518 double bin_width = 0.5;
519 double hist_low_limit = 0.0;
526 double theta_err = 1.0/57.3;
530 double tan_alpha_max = tan(theta_max + theta_err)/r0;
531 double tan_alpha_min = tan(theta_min - theta_err)/r0;
532 if(tan_alpha_min<0.0)tan_alpha_min=0.0;
533 for(
unsigned int i=0; i<seed.
hits.size(); i++){
535 if(!trkhit->
flags&VALID_HIT)
continue;
536 if(!trkhit->
flags&ON_CIRCLE)
continue;
537 if(!trkhit->
flags&IN_THETA_RANGE)
continue;
540 double q_sign = seed.
q>0.0 ? +1.0:-1.0;
541 double z_hit = trkhit->
hit->
pos.Z();
542 double zmin = z_hit - q_sign*trkhit->
phi_hit/tan_alpha_min;
543 double zmax = z_hit - q_sign*trkhit->
phi_hit/tan_alpha_max;
555 unsigned int imin = zmin<hist_low_limit ? 0:(
unsigned int)floor((zmin-hist_low_limit)/bin_width);
556 unsigned int imax = zmax<hist_low_limit ? 0:(
unsigned int)floor((zmax-hist_low_limit)/bin_width);
558 if(debug_level>3)
_DBG_<<
" -- phi_hit="<<trkhit->
phi_hit<<
" z_hit="<<z_hit<<endl;
559 if(debug_level>3)
_DBG_<<
" -- zmin="<<zmin<<
" zmax="<<zmax<<endl;
560 if(debug_level>3)
_DBG_<<
" -- imin="<<imin<<
" imax="<<imax<<endl;
564 if(imax<=0 || imin>=Nbins)
continue;
567 if(imin>=Nbins)imin=Nbins-1;
568 if(imax>=Nbins)imax=Nbins-1;
571 for(
unsigned int j=imin; j<=imax; j++)hist[j]++;
575 unsigned int istart=(
unsigned int)((TARGET_Z_MIN-hist_low_limit)/bin_width);
577 for(
unsigned int i=1; i<
Nbins; i++){
580 double z = hist_low_limit + bin_width*(0.5+(double)i);
581 if(z<TARGET_Z_MIN || z>TARGET_Z_MAX)
continue;
583 if(hist[i]>hist[istart]){
585 if(debug_level>3)
_DBG_<<
" -- istart="<<istart<<
" (z="<<hist_low_limit + bin_width*(0.5+(double)istart)<<
" , N="<<hist[i]<<
")"<<endl;
587 if(hist[i] == hist[istart])iend = i;
591 if(hist[istart]==0.0){
592 if(debug_level>3)
_DBG_<<
" - No entries in z-vertex hist. Seed aborted."<<endl;
597 seed.
z_min = hist_low_limit + bin_width*(0.5+(double)istart);
598 seed.
z_max = hist_low_limit + bin_width*(0.5+(double)iend);
600 if(debug_level>3)
_DBG_<<
"istart="<<istart<<
" iend="<<iend<<
" z_min="<<seed.
z_min<<
" z_max="<<seed.
z_max<<
" hits[istart]="<<hist[istart]<<endl;
604 for(
unsigned int i=0; i<seed.
hits.size(); i++){
606 trkhit->
flags &= ~IN_Z_RANGE;
607 if(!trkhit->
flags&VALID_HIT)
continue;
608 if(!trkhit->
flags&ON_CIRCLE)
continue;
609 if(!trkhit->
flags&IN_THETA_RANGE)
continue;
612 trkhit->
flags |= IN_Z_RANGE;
void setMomentum(const DVector3 &aMomentum)
const DFDCIntersection * hit
unsigned int NumAvailableHits(void)
virtual jerror_t init(void)
virtual jerror_t brun(JEventLoop *loop, int32_t runnumber)
void FindTheta(DFDCSeed &seed, double target_z_min, double target_z_max)
void FindSeeds(vector< DFDCSeed > &seeds)
void FindZ(DFDCSeed &seed, double theta_min, double theta_max)
vector< DFDCTrkHit * > hits
bool FDCSortByZincreasing(DTrackCandidate_factory_FDC::DFDCTrkHit *const &hit1, DTrackCandidate_factory_FDC::DFDCTrkHit *const &hit2)
double Dist2(const DFDCTrkHit *trkhit)
void FindThetaZ(DFDCSeed &seed)
virtual jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)
Invoked via JEventProcessor virtual method.
virtual jerror_t fini(void)
Invoked via JEventProcessor virtual method.
<A href="index.html#legend"> <IMG src="CORE.png" width="100"> </A>
void setPID(Particle_t locPID)
void setPosition(const DVector3 &aPosition)
void GetTrkHits(JEventLoop *loop)
void FillSeedHits(DFDCSeed &seed)
DTrackCandidate_factory_FDC()