Skip to content

Commit

Permalink
Merge pull request #1164 from JeffersonLab/TrackFinderFixes2
Browse files Browse the repository at this point in the history 
Track finder fixes2
  • Loading branch information
@sdobbs
sdobbs authored Jul 7, 2018
2 parents 9a726b3 + c6cfd06 commit cad2113
Show file tree
Hide file tree
Showing 18 changed files with 794 additions and 590 deletions.
84 changes: 63 additions & 21 deletions src/libraries/FDC/DFDCSegment_factory.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,9 +34,18 @@ using namespace std;
//}

bool DFDCSegment_package_cmp(const DFDCPseudo* a, const DFDCPseudo* b) {
return a->wire->layer>b->wire->layer;
if(a->wire->layer == b->wire->layer) {
if (fabs(a->time-b->time)<EPS){
double tsum_1=a->t_u+a->t_v;
double tsum_2=b->t_u+b->t_v;
return (tsum_1<tsum_2);
}
return(a->time<b->time);
}
return a->wire->layer>b->wire->layer;
}


DFDCSegment_factory::DFDCSegment_factory() {
_log = new JStreamLog(std::cout, "FDCSEGMENT >>");
*_log << "File initialized." << endMsg;
Expand Down Expand Up @@ -102,7 +111,7 @@ jerror_t DFDCSegment_factory::evnt(JEventLoop* eventLoop, uint64_t eventNo) {

// Find the segments in each package
for (int j=0;j<4;j++){
std::sort(package[j].begin(),package[j].end(),DFDCSegment_package_cmp);
std::stable_sort(package[j].begin(),package[j].end(),DFDCSegment_package_cmp);
// We need at least 3 points to define a circle, so skip if we don't
// have enough points.
if (package[j].size()>2) FindSegments(package[j]);
Expand All @@ -116,7 +125,7 @@ jerror_t DFDCSegment_factory::evnt(JEventLoop* eventLoop, uint64_t eventNo) {
// the dip angle and the z position of the closest approach to the beam line.
// Also returns predicted positions along the helical path.
//
jerror_t DFDCSegment_factory::RiemannLineFit(vector<const DFDCPseudo *>points,
jerror_t DFDCSegment_factory::RiemannLineFit(vector<const DFDCPseudo *>&points,
DMatrix &CR,vector<xyz_t>&XYZ){
unsigned int n=points.size()+1;
vector<int>bad(n); // Keep track of "bad" intersection points
Expand Down Expand Up @@ -161,10 +170,12 @@ jerror_t DFDCSegment_factory::RiemannLineFit(vector<const DFDCPseudo *>points,
double z=0,zlast=0;
double two_rc=2.*rc;
DVector2 oldxy=XYZ[0].xy;
unsigned int num_z=0;
for (unsigned int k=0;k<n;k++){
zlast=z;
z=XYZ[k].z;
if (fabs(z-zlast)<0.01) continue;
num_z++;

DVector2 diffxy=XYZ[k].xy-oldxy;
double var=CR(k,k);
Expand All @@ -188,7 +199,7 @@ jerror_t DFDCSegment_factory::RiemannLineFit(vector<const DFDCPseudo *>points,
//oldx=XYZ(k,0);
//oldy=XYZ(k,1);
oldxy=XYZ[k].xy;
}
}

Delta=sumv*sumxx-sumx*sumx;
double denom=sumv*sumxy-sumy*sumx;
Expand All @@ -215,6 +226,10 @@ jerror_t DFDCSegment_factory::RiemannLineFit(vector<const DFDCPseudo *>points,
else tanl=(zlast-zvertex)/sperp;
var_tanl=100.; // guess for now
}
if (num_z<3){
//_DBG_ << "Too few z points..." << endl;
return VALUE_OUT_OF_RANGE;
}
if (got_bad_intersection) return VALUE_OUT_OF_RANGE;
return NOERROR;
}
Expand Down Expand Up @@ -321,7 +336,7 @@ jerror_t DFDCSegment_factory::UpdatePositionsAndCovariance(unsigned int n,
// task of fitting points in (x,y) to a circle is transormed to the taks of
// fitting planes in (x,y, w=x^2+y^2) space
//
jerror_t DFDCSegment_factory::RiemannCircleFit(vector<const DFDCPseudo *>points,
jerror_t DFDCSegment_factory::RiemannCircleFit(vector<const DFDCPseudo *>&points,
DMatrix &CRPhi){
unsigned int n=points.size()+1;
DMatrix X(n,3);
Expand Down Expand Up @@ -438,11 +453,12 @@ jerror_t DFDCSegment_factory::RiemannCircleFit(vector<const DFDCPseudo *>points,
// length versus z. Uses RiemannCircleFit and RiemannLineFit.
//
jerror_t
DFDCSegment_factory::RiemannHelicalFit(vector<const DFDCPseudo*>points){
DFDCSegment_factory::RiemannHelicalFit(vector<const DFDCPseudo*>&points){
double Phi;
unsigned int num_measured=points.size();
unsigned int last_index=num_measured;
unsigned int num_points=num_measured+1;
unsigned int num_points=num_measured+1;
Ndof=num_points-3;
// list of points on track and corresponding covariance matrices
vector<xyz_t>XYZ(num_points);
DMatrix CR(num_points,num_points);
Expand Down Expand Up @@ -611,8 +627,7 @@ DFDCSegment_factory::RiemannHelicalFit(vector<const DFDCPseudo*>points){
rotation_sense=rotation_sense_;
rc=rc_,xc=xc_,yc=yc_,tanl=tanl_,zvertex=zvertex_,Phi1=Phi1_;
}
Ndof=num_points-3;


return NOERROR;
}

Expand All @@ -621,7 +636,7 @@ DFDCSegment_factory::RiemannHelicalFit(vector<const DFDCPseudo*>points){
// Associate nearest neighbors within a package with track segment candidates.
// Provide guess for (seed) track parameters
//
jerror_t DFDCSegment_factory::FindSegments(vector<const DFDCPseudo*>points){
jerror_t DFDCSegment_factory::FindSegments(vector<const DFDCPseudo*>&points){
// Create a vector to store whether or not a hit has been used
vector<bool>used(points.size());
unsigned int total_num_used=0;
Expand Down Expand Up @@ -716,7 +731,7 @@ jerror_t DFDCSegment_factory::FindSegments(vector<const DFDCPseudo*>points){
}
// Sort if we added another point
if (do_sort)
std::sort(neighbors.begin(),neighbors.end(),DFDCSegment_package_cmp);
std::stable_sort(neighbors.begin(),neighbors.end(),DFDCSegment_package_cmp);

// Arc lengths in helical model are referenced relative to the plane
// ref_plane within a segment. For a 6 hit segment, ref_plane=2 is
Expand All @@ -732,10 +747,13 @@ jerror_t DFDCSegment_factory::FindSegments(vector<const DFDCPseudo*>points){
// how well the points are defined at this stage. Use a simple circle
// fit assuming the circle goes through the origin.
if (rc<1.0) {
CircleFit(neighbors);
LineFit(neighbors);
if (CircleFit(neighbors)==NOERROR){
if (LineFit(neighbors)==NOERROR){
chisq=ComputeCircleChiSq(neighbors);
}
}
}

// Create a new segment
DFDCSegment *segment = new DFDCSegment;
segment->hits=neighbors;
Expand All @@ -746,8 +764,11 @@ jerror_t DFDCSegment_factory::FindSegments(vector<const DFDCPseudo*>points){
}
else {
// Fit assuming particle came from (x,y)=(0,0)
CircleFit(neighbors);
LineFit(neighbors);
if (CircleFit(neighbors)==NOERROR){
if (LineFit(neighbors)==NOERROR){
chisq=ComputeCircleChiSq(neighbors);
}
}

// Create a new segment
DFDCSegment *segment = new DFDCSegment;
Expand Down Expand Up @@ -796,17 +817,20 @@ jerror_t DFDCSegment_factory::FindSegments(vector<const DFDCPseudo*>points){
} // check if hit used already
}// loop over hits
if (added_hits){
sort(segment->hits.begin(),segment->hits.end(),DFDCSegment_package_cmp);
stable_sort(segment->hits.begin(),segment->hits.end(),DFDCSegment_package_cmp);
if (RiemannHelicalFit(segment->hits)==NOERROR){
// Since the cell size is 0.5 cm and the Lorentz deflection can be
// mm scale, a circle radius on the order of 1 cm is at the level of
// how well the points are defined at this stage. Use a simple circle
// fit assuming the circle goes through the origin.
if (rc<1.0) {
CircleFit(segment->hits);
LineFit(segment->hits);
if (CircleFit(segment->hits)==NOERROR){
if (LineFit(segment->hits)==NOERROR){
chisq=ComputeCircleChiSq(segment->hits);
}
}
}

FillSegmentData(segment);
}
}
Expand Down Expand Up @@ -941,11 +965,29 @@ jerror_t DFDCSegment_factory::CircleFit(vector<const DFDCPseudo *>&points){
if(fabs(denom)<1.0E-20)return UNRECOVERABLE_ERROR;
xc = (deltax*beta-deltay*gamma)/denom;
yc = (deltay*alpha-deltax*gamma)/denom;
rc = sqrt(xc*xc + yc*yc);
rc = sqrt(xc*xc + yc*yc);
Ndof=points.size()-2;

return NOERROR;
}

// Crude chi^2 calculation for circle fit that was forced to go through (0,0)
double DFDCSegment_factory::ComputeCircleChiSq(vector<const DFDCPseudo *>&neighbors) {
chisq=0.;
Phi1=atan2(neighbors[0]->xy.Y()-yc,neighbors[0]->xy.X()-xc);
double z0=neighbors[0]->wire->origin.z();
for (unsigned int m=0;m<neighbors.size();m++){
double sperp=rotation_sense*(neighbors[m]->wire->origin.z()-z0)/tanl;
double phi_s=Phi1+sperp/rc;
DVector2 XY(xc+rc*cos(phi_s),yc+rc*sin(phi_s));
chisq+=(XY-neighbors[m]->xy).Mod2();
// Very crude estimate: we assumed errors of size 1. in circle fit...
}

return chisq;
}


// Put fit results into the segment
void DFDCSegment_factory::FillSegmentData(DFDCSegment *segment){
// the particle's charge
Expand Down
9 changes: 5 additions & 4 deletions src/libraries/FDC/DFDCSegment_factory.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,14 +53,14 @@ class DFDCSegment_factory : public JFactory<DFDCSegment> {
}xyz_t;


jerror_t FindSegments(vector<const DFDCPseudo*>points);
jerror_t FindSegments(vector<const DFDCPseudo*>&points);
// jerror_t CorrectPoints(vector<DFDCPseudo*>point,DMatrix XYZ);
jerror_t GetHelicalTrackPosition(double z,const DFDCSegment *segment,
double &xpos,double &ypos);
jerror_t RiemannHelicalFit(vector<const DFDCPseudo*>points);
jerror_t RiemannCircleFit(vector<const DFDCPseudo*>points,
jerror_t RiemannHelicalFit(vector<const DFDCPseudo*>&points);
jerror_t RiemannCircleFit(vector<const DFDCPseudo*>&points,
DMatrix &CRPhi);
jerror_t RiemannLineFit(vector<const DFDCPseudo *>points,
jerror_t RiemannLineFit(vector<const DFDCPseudo *>&points,
DMatrix &CR,vector<xyz_t>&XYZ);
jerror_t UpdatePositionsAndCovariance(unsigned int n,double r1sq,
vector<xyz_t> &XYZ,DMatrix &CRPhi,
Expand All @@ -69,6 +69,7 @@ class DFDCSegment_factory : public JFactory<DFDCSegment> {
DMatrix &CRPhi, vector<const DFDCPseudo *>&points);
jerror_t CircleFit(vector<const DFDCPseudo *>&points);
jerror_t LineFit(vector<const DFDCPseudo *>&points);
double ComputeCircleChiSq(vector<const DFDCPseudo *>&neighbors);

void FillSegmentData(DFDCSegment *segment);

Expand Down
9 changes: 9 additions & 0 deletions src/libraries/TOF/DTOFPoint_factory.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,10 @@ bool Compare_TOFSpacetimeHitMatches_Distance(const DTOFPoint_factory::tof_spacet
return (locTOFSpacetimeHitMatch1.delta_r < locTOFSpacetimeHitMatch2.delta_r);
};

bool Compare_TOFPoint_Time(const DTOFPoint *locTOFPoint1, const DTOFPoint *locTOFPoint2) {
return locTOFPoint1->t < locTOFPoint2->t;
}

//------------------
// brun
//------------------
Expand Down Expand Up @@ -162,6 +166,11 @@ jerror_t DTOFPoint_factory::evnt(JEventLoop *loop, uint64_t eventnumber)
for(; locSetIterator != locUnusedTOFSpacetimeHits.end(); ++locSetIterator)
Create_UnMatchedTOFPoint(*locSetIterator);

// make sure all the hits are sorted by time (why not?)
// this helps with reproducibiliy problems...
std::sort(_data.begin(), _data.end(), Compare_TOFPoint_Time);


return NOERROR;
}

Expand Down
116 changes: 113 additions & 3 deletions src/libraries/TRACKING/DReferenceTrajectory.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,6 +314,118 @@ void DReferenceTrajectory::FastSwim(const DVector3 &pos, const DVector3 &mom,
}
}

// Faster version of swimmer that only deals with the region inside the bore of
// the magnet. This is intended to be used with the hit selector, so is a
// paired-down version of the usual trajectory.
void DReferenceTrajectory::FastSwimForHitSelection(const DVector3 &pos, const DVector3 &mom, double q){
DMagneticFieldStepper stepper(bfield, q, &pos, &mom);
if(step_size>0.0)stepper.SetStepSize(step_size);

// Step until we leave the active tracking region
swim_step_t *swim_step = this->swim_steps;
Nswim_steps = 0;
double itheta02 = 0.0;
double itheta02s = 0.0;
double itheta02s2 = 0.0;
double X0sum=0.0;
swim_step_t *last_step=NULL;

for(double s=0; fabs(s)<1000.; Nswim_steps++, swim_step++){

if(Nswim_steps>=this->max_swim_steps){
if (debug_level>0){
jerr<<__FILE__<<":"<<__LINE__<<" Too many steps in trajectory. Truncating..."<<endl;
}
break;
}

stepper.GetDirs(swim_step->sdir, swim_step->tdir, swim_step->udir);
stepper.GetPosMom(swim_step->origin, swim_step->mom);
swim_step->Ro = stepper.GetRo();
swim_step->s = s;
swim_step->t = 0.; // we do not need the time for our purpose...

// Magnetic field at current position
bfield->GetField(swim_step->origin,swim_step->B);

//magnitude of momentum and beta
double p_sq=swim_step->mom.Mag2();
double one_over_beta_sq=1.+mass_sq/p_sq;

// Add material if geom is not NULL
double dP = 0.0;
double dP_dx=0.0;
if(geom){
double KrhoZ_overA=0.0;
double rhoZ_overA=0.0;
double LogI=0.0;
double X0=0.0;
jerror_t err = geom->FindMatALT1(swim_step->origin, swim_step->mom, KrhoZ_overA, rhoZ_overA,LogI, X0);
if(err == NOERROR){
if(X0>0.0){
double p=sqrt(p_sq);
double delta_s = s;
if(last_step)delta_s -= last_step->s;
double radlen = delta_s/X0;

if(radlen>1.0E-5){ // PDG 2008 pg 271, second to last paragraph

// double theta0 = 0.0136*sqrt(one_over_beta_sq)/p*sqrt(radlen)*(1.0+0.038*log(radlen)); // From PDG 2008 eq 27.12
//double theta02 = theta0*theta0;
double factor=1.0+0.038*log(radlen);
double theta02=1.8496e-4*factor*factor*radlen*one_over_beta_sq/p_sq;

itheta02 += theta02;
itheta02s += s*theta02;
itheta02s2 += s*s*theta02;
X0sum+=X0;
}

// Calculate momentum loss due to ionization
dP_dx = dPdx(p, KrhoZ_overA, rhoZ_overA,LogI);
}
}
last_step = swim_step;
}
swim_step->itheta02 = itheta02;
swim_step->itheta02s = itheta02s;
swim_step->itheta02s2 = itheta02s2;
swim_step->invX0=Nswim_steps/X0sum;

if(step_size<0.0){ // step_size<0 indicates auto-calculated step size
// Adjust step size to take smaller steps in regions of high momentum loss
double my_step_size = 0.0001/fabs(dP_dx);
if(my_step_size>MAX_STEP_SIZE)my_step_size=MAX_STEP_SIZE; // maximum step size in cm
if(my_step_size<MIN_STEP_SIZE)my_step_size=MIN_STEP_SIZE; // minimum step size in cm

stepper.SetStepSize(my_step_size);
}

// Swim to next
double ds=stepper.FastStep(swim_step->B);

// Calculate momentum loss due to the step we're about to take
dP = ds*dP_dx;

// Adjust momentum due to ionization losses
if(dP!=0.0){
DVector3 pos, mom;
stepper.GetPosMom(pos, mom);
double ptot = mom.Mag() - dP; // correct for energy loss
if (ptot<0.005) {Nswim_steps++; break;}
mom.SetMag(ptot);
stepper.SetStartingParams(q, &pos, &mom);
}
s += ds;

double Rsq=swim_step->origin.Perp2();
double z=swim_step->origin.Z();
if (z>345. || z<17. || Rsq>60.*60.){
Nswim_steps++; break;
}
}
}

// Faster version of the swimmer that uses an alternate stepper and does not
// check for material boundaries.
void DReferenceTrajectory::FastSwim(const DVector3 &pos, const DVector3 &mom, double q,double smax, double zmin,double zmax){
Expand Down Expand Up @@ -472,6 +584,7 @@ void DReferenceTrajectory::FastSwim(const DVector3 &pos, const DVector3 &mom, do
hit_fcal=true;
}


// Exit the loop if we are already inside the volume of the BCAL
// and the radius is decreasing
if (Rsq<old_radius_sq && Rsq>Rsqmax_interior && z<407.0 && z>-100.0){
Expand Down Expand Up @@ -1401,9 +1514,6 @@ double DReferenceTrajectory::DistToRT(DVector3 hit, double *s,
return numeric_limits<double>::quiet_NaN();
}

if (std::isnan(Ro))
{
}
}

// Calculate distance along track ("s")
Expand Down
Loading

0 comments on commit cad2113

Please sign in to comment.