diff --git a/src/libraries/FDC/DFDCSegment_factory.cc b/src/libraries/FDC/DFDCSegment_factory.cc index 6f3708e3c..2996b9479 100644 --- a/src/libraries/FDC/DFDCSegment_factory.cc +++ b/src/libraries/FDC/DFDCSegment_factory.cc @@ -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)t_u+a->t_v; + double tsum_2=b->t_u+b->t_v; + return (tsum_1timetime); + } + return a->wire->layer>b->wire->layer; } + DFDCSegment_factory::DFDCSegment_factory() { _log = new JStreamLog(std::cout, "FDCSEGMENT >>"); *_log << "File initialized." << endMsg; @@ -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]); @@ -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(vectorpoints, +jerror_t DFDCSegment_factory::RiemannLineFit(vector&points, DMatrix &CR,vector&XYZ){ unsigned int n=points.size()+1; vectorbad(n); // Keep track of "bad" intersection points @@ -161,10 +170,12 @@ jerror_t DFDCSegment_factory::RiemannLineFit(vectorpoints, 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;kpoints, //oldx=XYZ(k,0); //oldy=XYZ(k,1); oldxy=XYZ[k].xy; - } + } Delta=sumv*sumxx-sumx*sumx; double denom=sumv*sumxy-sumy*sumx; @@ -215,6 +226,10 @@ jerror_t DFDCSegment_factory::RiemannLineFit(vectorpoints, 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; } @@ -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(vectorpoints, +jerror_t DFDCSegment_factory::RiemannCircleFit(vector&points, DMatrix &CRPhi){ unsigned int n=points.size()+1; DMatrix X(n,3); @@ -438,11 +453,12 @@ jerror_t DFDCSegment_factory::RiemannCircleFit(vectorpoints, // length versus z. Uses RiemannCircleFit and RiemannLineFit. // jerror_t -DFDCSegment_factory::RiemannHelicalFit(vectorpoints){ +DFDCSegment_factory::RiemannHelicalFit(vector&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 vectorXYZ(num_points); DMatrix CR(num_points,num_points); @@ -611,8 +627,7 @@ DFDCSegment_factory::RiemannHelicalFit(vectorpoints){ rotation_sense=rotation_sense_; rc=rc_,xc=xc_,yc=yc_,tanl=tanl_,zvertex=zvertex_,Phi1=Phi1_; } - Ndof=num_points-3; - + return NOERROR; } @@ -621,7 +636,7 @@ DFDCSegment_factory::RiemannHelicalFit(vectorpoints){ // Associate nearest neighbors within a package with track segment candidates. // Provide guess for (seed) track parameters // -jerror_t DFDCSegment_factory::FindSegments(vectorpoints){ +jerror_t DFDCSegment_factory::FindSegments(vector&points){ // Create a vector to store whether or not a hit has been used vectorused(points.size()); unsigned int total_num_used=0; @@ -716,7 +731,7 @@ jerror_t DFDCSegment_factory::FindSegments(vectorpoints){ } // 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 @@ -732,10 +747,13 @@ jerror_t DFDCSegment_factory::FindSegments(vectorpoints){ // 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; @@ -746,8 +764,11 @@ jerror_t DFDCSegment_factory::FindSegments(vectorpoints){ } 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; @@ -796,17 +817,20 @@ jerror_t DFDCSegment_factory::FindSegments(vectorpoints){ } // 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); } } @@ -941,11 +965,29 @@ jerror_t DFDCSegment_factory::CircleFit(vector&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&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;mwire->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 diff --git a/src/libraries/FDC/DFDCSegment_factory.h b/src/libraries/FDC/DFDCSegment_factory.h index 168aa71e2..beffc2e21 100644 --- a/src/libraries/FDC/DFDCSegment_factory.h +++ b/src/libraries/FDC/DFDCSegment_factory.h @@ -53,14 +53,14 @@ class DFDCSegment_factory : public JFactory { }xyz_t; - jerror_t FindSegments(vectorpoints); + jerror_t FindSegments(vector&points); // jerror_t CorrectPoints(vectorpoint,DMatrix XYZ); jerror_t GetHelicalTrackPosition(double z,const DFDCSegment *segment, double &xpos,double &ypos); - jerror_t RiemannHelicalFit(vectorpoints); - jerror_t RiemannCircleFit(vectorpoints, + jerror_t RiemannHelicalFit(vector&points); + jerror_t RiemannCircleFit(vector&points, DMatrix &CRPhi); - jerror_t RiemannLineFit(vectorpoints, + jerror_t RiemannLineFit(vector&points, DMatrix &CR,vector&XYZ); jerror_t UpdatePositionsAndCovariance(unsigned int n,double r1sq, vector &XYZ,DMatrix &CRPhi, @@ -69,6 +69,7 @@ class DFDCSegment_factory : public JFactory { DMatrix &CRPhi, vector&points); jerror_t CircleFit(vector&points); jerror_t LineFit(vector&points); + double ComputeCircleChiSq(vector&neighbors); void FillSegmentData(DFDCSegment *segment); diff --git a/src/libraries/TOF/DTOFPoint_factory.cc b/src/libraries/TOF/DTOFPoint_factory.cc index a85729178..55e5926ad 100644 --- a/src/libraries/TOF/DTOFPoint_factory.cc +++ b/src/libraries/TOF/DTOFPoint_factory.cc @@ -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 //------------------ @@ -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; } diff --git a/src/libraries/TRACKING/DReferenceTrajectory.cc b/src/libraries/TRACKING/DReferenceTrajectory.cc index 02d1b46e3..590fa171b 100644 --- a/src/libraries/TRACKING/DReferenceTrajectory.cc +++ b/src/libraries/TRACKING/DReferenceTrajectory.cc @@ -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..."<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_sizeB); + + // 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){ @@ -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 (RsqRsqmax_interior && z<407.0 && z>-100.0){ @@ -1401,9 +1514,6 @@ double DReferenceTrajectory::DistToRT(DVector3 hit, double *s, return numeric_limits::quiet_NaN(); } - if (std::isnan(Ro)) - { - } } // Calculate distance along track ("s") diff --git a/src/libraries/TRACKING/DReferenceTrajectory.h b/src/libraries/TRACKING/DReferenceTrajectory.h index 7189cfcc3..fca02c523 100644 --- a/src/libraries/TRACKING/DReferenceTrajectory.h +++ b/src/libraries/TRACKING/DReferenceTrajectory.h @@ -103,7 +103,7 @@ class DReferenceTrajectory{ swim_step_t* FindPlaneCrossing(const DVector3 &origin, DVector3 norm,int first_i=0, DetectorSystem_t detector=SYS_NULL) const; void Swim(const DVector3 &pos, const DVector3 &mom, double q=-1000.0,const TMatrixFSym *cov=NULL, double smax=2000.0, const DCoordinateSystem *wire=NULL); void FastSwim(const DVector3 &pos, const DVector3 &mom, double q,double smax=2000.0, double zmin=-100.,double zmax=1000.0); - + void FastSwimForHitSelection(const DVector3 &pos, const DVector3 &mom, double q); void FastSwim(const DVector3 &pos, const DVector3 &mom, DVector3 &last_pos, DVector3 &last_mom, diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index 833cbac94..42ce6e2ba 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory.cc @@ -23,6 +23,11 @@ /// Method 6: Tries to match stray unused cdc hits with fdc candidates /// Method 7: Alternate method to match leftover fdc candidates that /// have already been matched to other track candidates +/// Method 8: Attempts to "improve" a cdc candidate to be matched to an +/// FDC candidate with an assumption of the hit position in +/// outermost stereo straw +/// Methods 9-13: Various tricks to try to associate segments in the FDC +/// with other segments or tracks /// If a match is found, the code attempts to improve the track parameters by /// redoing the Riemann Circle fit with the additional hits. If an fdc /// candidate is matched to a cdc candidate, or several previously unused hits @@ -286,25 +291,32 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) const DTrackCandidate *srccan = cdctrackcandidates[i]; DVector3 mom=srccan->momentum(); DVector3 pos=srccan->position(); - double theta=mom.Theta(); - + // Propagate track to CDC endplate bool isForward=false; - if (theta0){ - // First do a quick projection using a helical model to see if it is - // worth adding this cdc candidate to the list of forward-going tracks - // that could pass into the FDC... - ProjectHelixToZ(cdc_endplate.z(),srccan->charge(),mom,pos); - if (pos.Perp()<48.5){ - // do an actual swim to the cdc endplate - mom=srccan->momentum(); - pos=srccan->position(); - stepper->SetCharge(srccan->charge()); - stepper->SwimToPlane(pos,mom,cdc_endplate,norm,NULL); + if (fdctrackcandidates.size()>0){ + // Check for candidates that appear to go backwards but are actually + //going forwards... + if (mom.Theta()>M_PI_2 && !sc_pos.empty()){ + TryToFlipDirection(schits,mom,pos); + } + if (mom.Theta()charge(),mom,pos); + if (pos.Perp()<48.5){ - cdc_endplate_projections.push_back(pos); - cdc_forward_ids.push_back(i); - isForward=true; + // do an actual swim to the cdc endplate + mom=srccan->momentum(); + pos=srccan->position(); + stepper->SetCharge(srccan->charge()); + stepper->SwimToPlane(pos,mom,cdc_endplate,norm,NULL); + if (pos.Perp()<48.5){ + cdc_endplate_projections.push_back(pos); + cdc_forward_ids.push_back(i); + isForward=true; + } } } } @@ -316,7 +328,7 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) // Variables for candidate number accounting int num_forward_cdc_cands_remaining=cdc_forward_ids.size(); int num_fdc_cands_remaining=fdctrackcandidates.size(); - + // Loop through the list of FDC candidates looking for matches between the // CDC and the FDC in the transition region. if (num_forward_cdc_cands_remaining>0){ @@ -329,7 +341,7 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) const DTrackCandidate *fdccan = fdctrackcandidates[i]; vectorsegments; fdccan->GetT(segments); - sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); + stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); if (segments[0]->package!=0) continue; // Flag for matching status @@ -420,15 +432,24 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) const DTrackCandidate *cdccan = cdctrackcandidates[cdc_forward_ids[j]]; vectorcdchits; cdccan->GetT(cdchits); - sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); + stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); // circle parameters can->rc=cdccan->rc; can->xc=cdccan->xc; can->yc=cdccan->yc; - can->setMomentum(cdccan->momentum()); - can->setPosition(cdccan->position()); + DVector3 mom=cdccan->momentum(); + DVector3 pos=cdccan->position(); + + // Check for candidates that appear to go backwards but are actually + //going forwards... + if (mom.Theta()>M_PI_2 && !sc_pos.empty()){ + TryToFlipDirection(schits,mom,pos); + } + + can->setMomentum(mom); + can->setPosition(pos); can->setPID(cdccan->PID()); for (unsigned int n=0;nmomentum(); DVector3 pos=cdccan->position(); - // Use information from the start counter to refine some candidates that - // seem to be pointing backwards but are probably pointing forwards. In - // these cases the circle projection intersects with a start counter - // paddle but the z-position at the radius r just outside the start - // counter barrel region is downstream of the nose, which does not make - // sense for a particle that is actually heading upstream... + // Check for candidates that appear to go backwards but are actually + //going forwards... if (mom.Theta()>M_PI_2 && !sc_pos.empty()){ - double zsc=sc_pos[0][1].z(); - if (pos.z()>zsc){ - unsigned int best_sc_sector_id=0; - double dphi_min=1000.; - double cand_phi=pos.Phi(); - for (unsigned int k=0;ksector-1; - double dphi=cand_phi-sc_pos[sector_id][0].Phi(); - if (dphi<-M_PI) dphi+=2.*M_PI; - if (dphi>M_PI) dphi-=2.*M_PI; - - if (fabs(dphi) +z - mom.SetMagThetaPhi(mom.Mag(),M_PI-mom.Theta(),mom.Phi()); - pos=sc_pos[best_sc_sector_id][1]; - // _DBG_ << " Event " << eventnumber << endl; - } - } - + TryToFlipDirection(schits,mom,pos); } - DTrackCandidate *can = new DTrackCandidate; can->setMomentum(mom); can->setPosition(pos); @@ -503,7 +496,7 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) // Get the cdc hits and add them to the candidate vectorcdchits; cdccan->GetT(cdchits); - sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); + stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); for (unsigned int n=0;nused_cdc_indexes[n]]=1; can->AddAssociatedObject(cdchits[n]); @@ -608,31 +601,36 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) } } - // We should be left with only single-segment fdc candidates: + // We should be left with only single-segment fdc candidates. We try to + // connect them together using alternate fitting techniques, either by + // forcing the circle to originate from the target or at the other extreme + // not including the fake point at the origin. if (num_fdc_cands_remaining){ - for (unsigned int j=0;jsegments; - srccan->GetT(segments); - - // redo circle fits for segments, forcing the circles to go through (0,0) - if (MatchMethod9(j,srccan,segments[0],fdctrackcandidates, - forward_matches)){ - if (DEBUG_LEVEL>0) _DBG_ << "Matched FDC segments using method #9" << endl; - num_fdc_cands_remaining-=2; - } - // Redo circle fit assuming track is almost straight - else if (MatchMethod10(j,srccan,segments[0],fdctrackcandidates, - forward_matches)){ - if (DEBUG_LEVEL>0) _DBG_ << "Matched FDC segments using method #10" << endl; - num_fdc_cands_remaining-=2; + bool matched_stray_segments=MatchStraySegments(forward_matches, + num_fdc_cands_remaining); + if (num_fdc_cands_remaining && matched_stray_segments){ + for (unsigned int i=0;imomentum().Mag()<0.5){ + MatchMethod12(can,forward_matches,num_fdc_cands_remaining); + } } - else { - // Not much we can do here -- add to the final list of candidates + } + } + if (num_fdc_cands_remaining){ + // Not much we can do here -- add to the final list of candidates + for (unsigned int j=0;jsegments; + srccan->GetT(segments); + const DFDCSegment *segment=segments[0]; + DTrackCandidate *can = new DTrackCandidate; // circle parameters can->rc=srccan->rc; @@ -641,11 +639,11 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) can->Ndof=srccan->Ndof; can->chisq=srccan->chisq; - can->setPID(srccan->PID()); + can->setPID(srccan->PID()); can->setMomentum(srccan->momentum()); can->setPosition(srccan->position()); - for (unsigned int n=0;nhits.size();n++){ - const DFDCPseudo *fdchit=segments[0]->hits[n]; + for (unsigned int n=0;nhits.size();n++){ + const DFDCPseudo *fdchit=segment->hits[n]; can->AddAssociatedObject(fdchit); } @@ -673,7 +671,7 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) vectorusedfdchits; candidate->GetT(usedfdchits); // Sort the hits - sort(usedfdchits.begin(),usedfdchits.end(),FDCHitSortByLayerincreasing); + stable_sort(usedfdchits.begin(),usedfdchits.end(),FDCHitSortByLayerincreasing); if ((usedfdchits[0]->wire->layer-1)/6==0){ MatchMethod6(candidate,usedfdchits,used_cdc_hits,num_unmatched_cdcs); } @@ -691,7 +689,7 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) vectormycdchits; candidate->GetT(mycdchits); - if (mycdchits.size()>=6 || myfdchits.size()>=3){ + if (mycdchits.size()>=3 || myfdchits.size()>=3){ _data.push_back(candidate); } else delete candidate; @@ -1202,7 +1200,7 @@ bool DTrackCandidate_factory::MatchMethod1(const DTrackCandidate *fdccan, // JANA does not maintain the order that the segments were added // as associated objects. Therefore, we need to reorder them here // so segment[0] is the most upstream segment. - sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); + stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); unsigned int cdc_index=cdc_forward_ids[jmin]; const DTrackCandidate *cdccan=cdctrackcandidates[cdc_index]; @@ -1212,7 +1210,7 @@ bool DTrackCandidate_factory::MatchMethod1(const DTrackCandidate *fdccan, cdccan->GetT(cdchits); // Sort CDC hits by layer - sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); + stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); // Abort if the track would have to pass from one quadrant into another // quadrant (this is more likely two roughly back-to-back tracks rather than @@ -1221,7 +1219,7 @@ bool DTrackCandidate_factory::MatchMethod1(const DTrackCandidate *fdccan, const DVector3 cdc_wire_origin=cdchits[0]->wire->origin; if (cdc_wire_origin.x()*fdc_hit_pos.X()<0. || cdc_wire_origin.y()*fdc_hit_pos.Y()<0.){ - if (DEBUG_LEVEL>1) _DBG_ << "Skipping match of potential back-to-back tracks." <0) _DBG_ << "Skipping match of potential back-to-back tracks." <cdchits; cdccan->GetT(cdchits); - sort(cdchits.begin(),cdchits.end(),CDCHitSortByLayerincreasing); + stable_sort(cdchits.begin(),cdchits.end(),CDCHitSortByLayerincreasing); // Variables to count the number of matching hits and the match fraction unsigned int num_match=0; @@ -1357,7 +1355,7 @@ bool DTrackCandidate_factory::MatchMethod1(const DTrackCandidate *fdccan, // JANA does not maintain the order that the segments were added // as associated objects. Therefore, we need to reorder them here // so segment[0] is the most upstream segment. - sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); + stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); // Abort if the track would have to pass from one quadrant into the opposite // quadrant (this is more likely two roughly back-to-back tracks rather than @@ -1366,7 +1364,7 @@ bool DTrackCandidate_factory::MatchMethod1(const DTrackCandidate *fdccan, const DVector3 cdc_wire_origin=cdchits[0]->wire->origin; if (cdc_wire_origin.x()*fdc_hit_pos.X()<0. && cdc_wire_origin.y()*fdc_hit_pos.Y()<0.){ - if (DEBUG_LEVEL>1) _DBG_ << "Skipping match of potential back-to-back tracks." <0) _DBG_ << "Skipping match of potential back-to-back tracks." <cdchits; cdccan->GetT(cdchits); - sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); + stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); // loop over fdc candidates for (unsigned int k=0;ksegments; fdccan->GetT(segments); - sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); + stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); // Only try to match candidate if this fdc candidate has hits in the // first package @@ -1498,7 +1496,7 @@ bool DTrackCandidate_factory::MatchMethod1(const DTrackCandidate *fdccan, const DVector3 cdc_wire_origin=cdchits[0]->wire->origin; if (cdc_wire_origin.x()*fdc_hit_pos.X()<0. && cdc_wire_origin.y()*fdc_hit_pos.Y()<0.){ - if (DEBUG_LEVEL>1) _DBG_ << "Skipping match of potential back-to-back tracks." <0) _DBG_ << "Skipping match of potential back-to-back tracks." <GetT(src_segments); if (src_segments.size()==1) return false; - sort(src_segments.begin(), src_segments.end(), SegmentSortByLayerincreasing); + stable_sort(src_segments.begin(), src_segments.end(), SegmentSortByLayerincreasing); if (DEBUG_LEVEL>0) _DBG_ << "Attempting matching method #4..." <segments; fdccan->GetT(segments); - sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); + stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); const DFDCPseudo *firsthit=segments[0]->hits[0]; unsigned int pack2_first=segments[0]->package; @@ -1838,7 +1836,7 @@ bool DTrackCandidate_factory::MatchMethod5(DTrackCandidate *can, // JANA does not maintain the order that the segments were added // as associated objects. Therefore, we need to reorder them here // so segment[0] is the most upstream segment. - sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); + stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); // Abort if the track would have to pass from one quadrant into the opposite // quadrant (this is more likely two roughly back-to-back tracks rather than @@ -1847,7 +1845,7 @@ bool DTrackCandidate_factory::MatchMethod5(DTrackCandidate *can, const DVector3 cdc_wire_origin=cdchits[0]->wire->origin; if (cdc_wire_origin.x()*fdc_hit_pos.X()<0. && cdc_wire_origin.y()*fdc_hit_pos.Y()<0.){ - if (DEBUG_LEVEL>1) _DBG_ << "Skipping match of potential back-to-back tracks." <0) _DBG_ << "Skipping match of potential back-to-back tracks." <momentum(); double q=can->charge(); - // Variables to keep track of cdc hits - unsigned int num_match_cdc=0; - unsigned int num_axial=0; + // Variables to keep track of cdc hits + bool got_inner_index=false; + unsigned int inner_index=0; int id_for_smallest_dr=-1; int old_ring=-1; - vectormatched_ids; double dr2_min=1e6; for (unsigned int k=0;kwire->ring!=old_ring){ if (id_for_smallest_dr>=0){ if (!used_cdc_hits[id_for_smallest_dr]){ - num_match_cdc++; num_unmatched_cdcs--; used_cdc_hits[id_for_smallest_dr]=1; can->used_cdc_indexes.push_back(id_for_smallest_dr); can->AddAssociatedObject(mycdchits[id_for_smallest_dr]); - if (mycdchits[id_for_smallest_dr]->is_stereo==false){ - num_axial++; - matched_ids.push_back(id_for_smallest_dr); + if (got_inner_index==false){ + inner_index=id_for_smallest_dr; + got_inner_index=true; } } } @@ -1969,26 +1965,17 @@ void DTrackCandidate_factory::MatchMethod6(DTrackCandidate *can, // Grab the last potential match if (id_for_smallest_dr>=0){ if (!used_cdc_hits[id_for_smallest_dr]){ - num_match_cdc++; num_unmatched_cdcs--; used_cdc_hits[id_for_smallest_dr]=1; can->used_cdc_indexes.push_back(id_for_smallest_dr); can->AddAssociatedObject(mycdchits[id_for_smallest_dr]); - - if (mycdchits[id_for_smallest_dr]->is_stereo==false){ - num_axial++; - matched_ids.push_back(id_for_smallest_dr); - } } } // We found some matched hits. Update the start position of the candidate. - if (num_match_cdc>0){ + if (got_inner_index){ const DFDCPseudo *firsthit=fdchits[0]; - - int axial_id=-1; - if (num_axial>0) axial_id=matched_ids[0]; // Compute the average Bz double Bz_avg=0.,denom=0.; @@ -2012,13 +1999,21 @@ void DTrackCandidate_factory::MatchMethod6(DTrackCandidate *can, fit.y0=pos.y()+fit.h*fit.r0*cos(fit.phi); fit.tanl=tan(M_PI_2-mom.Theta()); fit.z_vertex=0; // this will be changed later + if (GetPositionAndMomentum(fit,Bz_avg,mycdchits[inner_index]->wire->origin,pos, + mom)==NOERROR){ + // Find the z-position at the new position in x and y + DVector2 xy0(pos.X(),pos.Y()); + double tworc=2.*fit.r0; + double ratio=(firsthit->xy-xy0).Mod()/tworc; + double sperp=(ratio<1.)?tworc*asin(ratio):tworc*M_PI_2; + pos.SetZ(firsthit->wire->origin.z()-sperp*fit.tanl); - // Update the position and momentum of the candidate based on the most - // recent fit results - UpdatePositionAndMomentum(can,firsthit,fit,Bz_avg,axial_id); + // update the track parameters + can->setMomentum(mom); + can->setPosition(pos); - if (DEBUG_LEVEL>0) _DBG_ << "... matched stray CDC hits ..." << endl; - + if (DEBUG_LEVEL>0) _DBG_ << "... matched stray CDC hits ..." << endl; + } } // match at least one cdc hit } @@ -2034,7 +2029,7 @@ bool DTrackCandidate_factory::MatchMethod7(DTrackCandidate *srccan, srccan->GetT(fdchits); if (fdchits.size()==0) return false; - sort(fdchits.begin(),fdchits.end(),FDCHitSortByLayerincreasing); + stable_sort(fdchits.begin(),fdchits.end(),FDCHitSortByLayerincreasing); unsigned int pack1_last=(fdchits[fdchits.size()-1]->wire->layer-1)/6; for (unsigned int k=0;ksegments; fdccan->GetT(segments); - sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); + stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); const DFDCPseudo *firsthit=segments[0]->hits[0]; unsigned int pack2_first=segments[0]->package; @@ -2102,6 +2097,7 @@ bool DTrackCandidate_factory::MatchMethod7(DTrackCandidate *srccan, // axial straws to the list of hits to use in the circle fit vectorcdchits; srccan->GetT(cdchits); + stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); for (unsigned int i=0;iis_stereo==false){ double cov=0.213; //guess @@ -2203,7 +2199,7 @@ bool DTrackCandidate_factory::MatchMethod8(const DTrackCandidate *cdccan, // Get the cdc hits associated with this track candidate vectorcdchits; cdccan->GetT(cdchits); - sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); + stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); unsigned int last_index=cdchits.size()-1; if (cdchits[last_index]->is_stereo==true @@ -2235,7 +2231,7 @@ bool DTrackCandidate_factory::MatchMethod8(const DTrackCandidate *cdccan, // Fit the points to a circle assuming that the circle goes through the // origin. - fit.FitCircle(); + if (fit.FitCircle()!=NOERROR) return false; // Determine the tangent of the dip angle double tworc=2.*fit.r0; @@ -2273,7 +2269,7 @@ bool DTrackCandidate_factory::MatchMethod8(const DTrackCandidate *cdccan, // Get the segment data vectorsegments; fdccan->GetT(segments); - sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); + stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); // only do this for candidates with segments in the first package if (segments[0]->package>0) continue; @@ -2434,7 +2430,7 @@ bool DTrackCandidate_factory::MatchMethod9(unsigned int src_index, const DFDCPseudo *hit=segment->hits[n]; fit1.AddHit(hit); } - fit1.FitCircle(); + if (fit1.FitCircle()!=NOERROR) return false; // Sense of rotation fit1.h=q*FactorForSenseOfRotation; @@ -2471,7 +2467,7 @@ bool DTrackCandidate_factory::MatchMethod9(unsigned int src_index, const DFDCPseudo *hit=segments2[0]->hits[n]; fit2.AddHit(hit); } - fit2.FitCircle(); + if (fit2.FitCircle()!=NOERROR) continue; // Match using centers of circles double dx=fit1.x0-fit2.x0; @@ -2569,7 +2565,7 @@ bool DTrackCandidate_factory::MatchMethod9(unsigned int src_index, trackcandidates.push_back(can); - if (DEBUG_LEVEL>1){ + if (DEBUG_LEVEL>0){ _DBG_ << "Match method #9 succeeded..." << endl; } @@ -2632,8 +2628,7 @@ bool DTrackCandidate_factory::MatchMethod10(unsigned int src_index, int pack2=segments2[0]->package; if (abs(pack1-pack2)>0){ - // Get hits from the second segment and redo fit forcing circle - // to go through (0,0) + // Get hits from the second segment and redo fit DHelicalFit fit2; for (unsigned int n=0;nhits.size();n++){ const DFDCPseudo *hit=segments2[0]->hits[n]; @@ -2726,7 +2721,7 @@ bool DTrackCandidate_factory::MatchMethod10(unsigned int src_index, trackcandidates.push_back(can); - if (DEBUG_LEVEL>1){ + if (DEBUG_LEVEL>0){ _DBG_ << "Match method #10 succeeded..." << endl; } return true; @@ -2888,20 +2883,20 @@ bool DTrackCandidate_factory::MatchMethod12(DTrackCandidate *can, // Get the segments belonging to the fdc candidate vectorsegments; fdccan->GetT(segments); - sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); + stable_sort(segments.begin(), segments.end(), SegmentSortByLayerincreasing); // Get the hits from the input candidate vectormyfdchits; can->GetT(myfdchits); if (myfdchits.size()>0){ - sort(myfdchits.begin(),myfdchits.end(),FDCHitSortByLayerincreasing); + stable_sort(myfdchits.begin(),myfdchits.end(),FDCHitSortByLayerincreasing); unsigned int last_package =(myfdchits[myfdchits.size()-1]->wire->layer-1)/6; // look for something downstream... if (last_packagepackage){ vectorcdchits; can->GetT(cdchits); - sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); + stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); // Variables for computing average Bz double Bz=0; unsigned int num_hits=0; @@ -3002,7 +2997,7 @@ bool DTrackCandidate_factory::MatchMethod12(DTrackCandidate *can, // Get the cdc hits belonging to the track vectorcdchits; can->GetT(cdchits); - sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); + stable_sort(cdchits.begin(), cdchits.end(), CDCHitSortByLayerincreasing); // Variables for computing average Bz double Bz=0; @@ -3093,6 +3088,159 @@ bool DTrackCandidate_factory::MatchMethod12(DTrackCandidate *can, return false; } +// Attempts to match two single-segment track candidates by only using the +// hit information to fit circles (i.e., not requiring the tracks to originate +// from the target region.) The idea is to recover some tracks arising from +// secondary interactions or decays beyond the target. +bool DTrackCandidate_factory::MatchMethod13(unsigned int src_index, + const DTrackCandidate *srccan, + const DFDCSegment *segment, + vector&cands, + vector &forward_matches){ + if (DEBUG_LEVEL>0){ + _DBG_ << "Attempting Match method #13..." << endl; + } + double q=srccan->charge(); + int pack1=segment->package; + + // Get hits from segment and redo fit + DHelicalFit fit1; + for (unsigned int n=0;nhits.size();n++){ + const DFDCPseudo *hit=segment->hits[n]; + fit1.AddHit(hit); + } + + if (fit1.FitCircleRiemann(srccan->rc)!=NOERROR) return false; + + // Guess for theta and z from input candidates + double theta=srccan->momentum().Theta(); + fit1.tanl=tan(M_PI_2-theta); + fit1.z_vertex=srccan->position().Z(); + + // Redo line fit + if (fit1.FitLineRiemann()!=NOERROR) return false; + + // Find the magnetic field at the first hit in the first segment + double x=segment->hits[0]->xy.X(); + double y=segment->hits[0]->xy.Y(); + double z=segment->hits[0]->wire->origin.z(); + double Bz=fabs(bfield->GetBz(x,y,z)); + + // Get position and momentum for the first segment + DVector3 mypos,mymom; + GetPositionAndMomentum(z,fit1,Bz,mypos,mymom); + + // Loop over the rest of the fdc track candidates, skipping those that have + // already been used + for (unsigned int k=src_index+1;ksegments2; + can2->GetT(segments2); + + int pack2=segments2[0]->package; + if (abs(pack1-pack2)>0){ + // Get hits from the second segment and redo fit + DHelicalFit fit2; + for (unsigned int n=0;nhits.size();n++){ + const DFDCPseudo *hit=segments2[0]->hits[n]; + fit2.AddHit(hit); + } + + if (fit2.FitCircleRiemann(segments2[0]->rc)!=NOERROR) continue; + + // Guess for theta and z from input candidates + theta=can2->momentum().Theta(); + fit2.tanl=tan(M_PI_2-theta); + fit2.z_vertex=srccan->position().Z(); + + // Redo line fit + if (fit2.FitLineRiemann()!=NOERROR) continue; + + // Try to match segments by swimming through the field + if (MatchMethod11(q,mypos,mymom,fit2,segment,segments2[0])){ + forward_matches[k]=1; + forward_matches[src_index]=1; + + // Create a new DTrackCandidate for output + DTrackCandidate *can = new DTrackCandidate; + + // variables for finding + double Bz=0; + int num_hits=0; + + // Add hits from first segment as associated objects + for (unsigned int n=0;nhits.size();n++){ + const DFDCPseudo *fdchit=segment->hits[n]; + can->AddAssociatedObject(fdchit); + + Bz+=bfield->GetBz(fdchit->xy.X(),fdchit->xy.Y(), + fdchit->wire->origin.z()); + num_hits++; + } + + // Add the hits from the second segment to the first fit object + // and refit the circle. Also add them as associated objects + // to the new candidate. + for (unsigned int n=0;nhits.size();n++){ + const DFDCPseudo *hit=segments2[0]->hits[n]; + fit1.AddHit(hit); + can->AddAssociatedObject(hit); + + Bz+=bfield->GetBz(hit->xy.X(),hit->xy.Y(), + hit->wire->origin.z()); + num_hits++; + } + Bz=fabs(Bz)/double(num_hits); + + // Initialize variables needed for output + DVector3 mom=srccan->momentum(); + DVector3 pos=srccan->position(); + double q=srccan->charge(); + can->Ndof=srccan->Ndof; + can->chisq=srccan->chisq; + + if (fit1.FitCircleRiemann(fit1.r0)==NOERROR){ + can->Ndof=fit1.ndof; + can->chisq=fit1.chisq; + + // Redo line fit + fit1.FitLineRiemann(); + + // Guess charge from fit + fit1.h=GetSenseOfRotation(fit1,segments2[0]->hits[0], + srccan->position()); + q=FactorForSenseOfRotation*fit1.h; + + // put z position just upstream of the first hit in z + const DHFHit_t *myhit=fit1.GetHits()[0]; + pos.SetXYZ(myhit->x,myhit->y,myhit->z); + GetPositionAndMomentum(myhit->z-1.,fit1,Bz,pos,mom); + } + // circle parameters + can->rc=fit1.r0; + can->xc=fit1.x0; + can->yc=fit1.y0; + + can->setPID((q > 0.0) ? PiPlus : PiMinus); + can->setPosition(pos); + can->setMomentum(mom); + + trackcandidates.push_back(can); + + if (DEBUG_LEVEL>0){ + _DBG_ << "Match method #13 succeeded..." << endl; + } + return true; + } // minimum number of matching hits + } // different packages? + } // already matched? + } // loop over tracks + + return false; +} + // Update the momentum and position entries for the candidate based on the // results of the most recent helical fit void DTrackCandidate_factory::UpdatePositionAndMomentum(DTrackCandidate *can, @@ -3174,3 +3322,93 @@ DTrackCandidate_factory::GetPositionAndMomentum(double z,DHelicalFit &fit, return NOERROR; } + +// Use information from the start counter to try to correct the momentum and +// position for tracks seem to be pointing backwards but are probably pointing +// forwards. In these cases the circle projection intersects with a start +// counter paddle but the z-position at the radius r just outside the start +// counter barrel region is downstream of the nose, which does not make +// sense for a particle that is actually heading upstream... +bool DTrackCandidate_factory::TryToFlipDirection(vector&schits, + DVector3 &mom,DVector3 &pos) const{ + if (schits.size()==0) return false; + if (DEBUG_LEVEL>0){ + _DBG_ << "Attempting to flip direction of track..." << endl; + } + + double zsc=sc_pos[0][1].z(); + if (pos.z()>zsc){ + unsigned int best_sc_sector_id=0; + double dphi_min=1000.; + double cand_phi=pos.Phi(); + for (unsigned int k=0;ksector-1; + double dphi=cand_phi-sc_pos[sector_id][0].Phi(); + if (dphi<-M_PI) dphi+=2.*M_PI; + if (dphi>M_PI) dphi-=2.*M_PI; + + if (fabs(dphi) +z + mom.SetMagThetaPhi(mom.Mag(),M_PI-mom.Theta(),mom.Phi()); + pos=sc_pos[best_sc_sector_id][1]; + // _DBG_ << " Event " << eventnumber << endl; + if (DEBUG_LEVEL>0){ + _DBG_ << "Changing direction of CDC candidate..." << endl; + } + return true; + } + } + + + return false; +} + +// Last ditch attempt to match stray segments in the FDC to other track stubs +// after all other matching techniques have been tried... +bool DTrackCandidate_factory::MatchStraySegments(vector &forward_matches, + int &num_fdc_cands_remaining){ + if (DEBUG_LEVEL>0){ + _DBG_ << "Attempting to match stray FDC segments..." << endl; + } + + bool got_new_match=false; + for (unsigned int j=0;jsegments; + srccan->GetT(segments); + const DFDCSegment *segment=segments[0]; + + // redo circle fits for segments, forcing the circles to go through (0,0) + if (MatchMethod9(j,srccan,segment,fdctrackcandidates,forward_matches)){ + if (DEBUG_LEVEL>0) _DBG_ << "Matched FDC segments using method #9" << endl; + num_fdc_cands_remaining-=2; + got_new_match=true; + } + // Redo circle fit assuming track is almost straight + else if (MatchMethod10(j,srccan,segment,fdctrackcandidates, + forward_matches)){ + if (DEBUG_LEVEL>0) _DBG_ << "Matched FDC segments using method #10" << endl; + num_fdc_cands_remaining-=2; + got_new_match=true; + } + // Redo circle fits without any assumption that the tracks originate + // from the target + else if (MatchMethod13(j,srccan,segment,fdctrackcandidates, + forward_matches)){ + if (DEBUG_LEVEL>0) _DBG_ << "Matched FDC segments using method #13" << endl; + num_fdc_cands_remaining-=2; + got_new_match=true; + } + } // not already matched + } + return got_new_match; +} diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.h b/src/libraries/TRACKING/DTrackCandidate_factory.h index f46cd188b..972ea4465 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.h +++ b/src/libraries/TRACKING/DTrackCandidate_factory.h @@ -128,8 +128,15 @@ class DTrackCandidate_factory:public JFactory{ DHelicalFit &fit2,const DFDCSegment *segment1, const DFDCSegment *segment2); bool MatchMethod12(DTrackCandidate *srccan,vector &forward_matches, - int &num_fdc_cands_remaining); - + int &num_fdc_cands_remaining); + bool MatchMethod13(unsigned int src_index,const DTrackCandidate *srccan, + const DFDCSegment *segment, + vector&cands, + vector &forward_matches); + bool TryToFlipDirection(vector&scihits, + DVector3 &mom,DVector3 &pos) const; + bool MatchStraySegments(vector &forward_matches, + int &num_fdc_cands_remaining); private: const DMagneticFieldMap *bfield; diff --git a/src/libraries/TRACKING/DTrackCandidate_factory_CDC.cc b/src/libraries/TRACKING/DTrackCandidate_factory_CDC.cc index b54afb6cc..b33900b02 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_CDC.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory_CDC.cc @@ -274,7 +274,7 @@ jerror_t DTrackCandidate_factory_CDC::evnt(JEventLoop *locEventLoop, uint64_t ev cout << "final fit track circles" << endl; Print_TrackCircles(locCDCTrackCircles); } - sort(locCDCTrackCircles.begin(), locCDCTrackCircles.end(), CDCSortByChiSqPerNDFDecreasing); //sort by circle-fit weighted chisq/ndf (largest first) + stable_sort(locCDCTrackCircles.begin(), locCDCTrackCircles.end(), CDCSortByChiSqPerNDFDecreasing); //sort by circle-fit weighted chisq/ndf (largest first) Handle_StereoAndFilter(locCDCTrackCircles, true); //true: final pass: don't need to filter any more, get improved theta-z (although will reject if bad/no theta/z) if(DEBUG_LEVEL > 5) @@ -462,7 +462,7 @@ jerror_t DTrackCandidate_factory_CDC::Get_CDCHits(JEventLoop* loop) // Sort the individual superlayer lists by decreasing values of R for(size_t i = 0; i < cdchits_by_superlayer.size(); ++i) - sort(cdchits_by_superlayer[i].begin(), cdchits_by_superlayer[i].end(), CDCSortByRdecreasing); + stable_sort(cdchits_by_superlayer[i].begin(), cdchits_by_superlayer[i].end(), CDCSortByRdecreasing); // Filter out noise hits. All hits are initially flagged as "noise". // Hits with a neighbor within MAX_HIT_DIST have their noise flags cleared. @@ -1622,7 +1622,7 @@ bool DTrackCandidate_factory_CDC::Build_TrackCircles(vector& l Fit_Circles(locCDCTrackCircles, false, false); if(locCDCTrackCircles.empty()) return true; - sort(locCDCTrackCircles.begin(), locCDCTrackCircles.end(), CDCSortByChiSqPerNDFDecreasing); //sort by circle-fit weighted chisq/ndf (largest first) + stable_sort(locCDCTrackCircles.begin(), locCDCTrackCircles.end(), CDCSortByChiSqPerNDFDecreasing); //sort by circle-fit weighted chisq/ndf (largest first) if(DEBUG_LEVEL > 5) { @@ -2848,7 +2848,7 @@ void DTrackCandidate_factory_CDC::Truncate_TrackCircles(vector //now fit the newly-truncated track circles Fit_Circles(locCDCTrackCircles, true, false); //true: fit only truncated circles //false: don't add stereo intersections - sort(locCDCTrackCircles.begin(), locCDCTrackCircles.end(), CDCSortByChiSqPerNDFDecreasing); //sort by fit chisq/ndf + stable_sort(locCDCTrackCircles.begin(), locCDCTrackCircles.end(), CDCSortByChiSqPerNDFDecreasing); //sort by fit chisq/ndf if(DEBUG_LEVEL > 5) { cout << "Post-SL7-turncation track circles" << endl; @@ -3023,7 +3023,7 @@ void DTrackCandidate_factory_CDC::Truncate_TrackCircles(vector //now fit the newly-truncated track circles Fit_Circles(locCDCTrackCircles, true, false); //true: fit only truncated circles //false: don't add stereo intersections - sort(locCDCTrackCircles.begin(), locCDCTrackCircles.end(), CDCSortByChiSqPerNDFDecreasing); //sort by fit chisq/ndf + stable_sort(locCDCTrackCircles.begin(), locCDCTrackCircles.end(), CDCSortByChiSqPerNDFDecreasing); //sort by fit chisq/ndf if(DEBUG_LEVEL > 5) { cout << "Post-SL4-turncation track circles" << endl; @@ -3527,7 +3527,7 @@ void DTrackCandidate_factory_CDC::Select_ThetaZStereoHits(const DCDCTrackCircle* unsigned int locSuperLayer = locSuperLayerSeeds[loc_k]->dSuperLayer; Calc_StereoHitDeltaPhis(locSuperLayer, locHitsBySuperLayer[locSuperLayer], locCDCTrackCircle, locDeltaPhis); } - sort(locDeltaPhis.begin(), locDeltaPhis.end(), CDCSort_DeltaPhis); + stable_sort(locDeltaPhis.begin(), locDeltaPhis.end(), CDCSort_DeltaPhis); if(locDeltaPhis.size() <= MIN_PRUNED_STEREO_HITS) { @@ -3849,7 +3849,7 @@ bool DTrackCandidate_factory_CDC::Find_ThetaZ_Regression(const DHelicalFit* locF } // Now, sort the entries - sort(intersections.begin(), intersections.end(), CDCSort_Intersections); + stable_sort(intersections.begin(), intersections.end(), CDCSort_Intersections); // Compute the arc lengths between the origin in x and y and (xi,yi) vector arclengths(intersections.size()); @@ -4258,7 +4258,7 @@ void DTrackCandidate_factory_CDC::Filter_TrackCircles_Stereo(vector 5) { @@ -4308,11 +4308,11 @@ void DTrackCandidate_factory_CDC::Filter_TrackCircles_Stereo(vector 5) { diff --git a/src/libraries/TRACKING/DTrackCandidate_factory_FDC.cc b/src/libraries/TRACKING/DTrackCandidate_factory_FDC.cc index eab72123b..f3a872847 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_FDC.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory_FDC.cc @@ -286,7 +286,7 @@ void DTrackCandidate_factory_FDC::FillSeedHits(DFDCSeed &seed) } // Sort hits by increasing z - sort(seed.hits.begin(), seed.hits.end(), FDCSortByZincreasing); + stable_sort(seed.hits.begin(), seed.hits.end(), FDCSortByZincreasing); // Below, we fill in the phi_hit field of the DFDCTrkHit objects for // this candidate. We do so incrementally to try and accomodate diff --git a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc index 6db0de27e..b8c06b428 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc @@ -18,11 +18,6 @@ #include #include -#define MAX_SEGMENTS 20 -#define HALF_PACKAGE 6.0 -#define FDC_OUTER_RADIUS 50.0 -#define BEAM_VAR 0.0208 // cm^2 -#define HIT_CHI2_CUT 10.0 /// /// DTrackCandidate_factory_FDCCathodes::brun(): /// @@ -39,18 +34,6 @@ jerror_t DTrackCandidate_factory_FDCCathodes::brun(JEventLoop* eventLoop, _DBG_<< "FDC geometry not available!" <GetCDCEndplate(endplate_z,endplate_dz,endplate_rmin,endplate_rmax); - endplate_z+=endplate_dz; - JCalibration *jcalib = dapp->GetJCalibration(runnumber); map targetparms; @@ -64,22 +47,8 @@ jerror_t DTrackCandidate_factory_FDCCathodes::brun(JEventLoop* eventLoop, DEBUG_HISTS=false; gPARMS->SetDefaultParameter("TRKFIND:DEBUG_HISTS", DEBUG_HISTS); - APPLY_MOMENTUM_CORRECTION=false; - gPARMS->SetDefaultParameter("TRKFIND:APPLY_MOMENTUM_CORRECTION",APPLY_MOMENTUM_CORRECTION); - p_factor1=1.61*M_PI/180.; - p_factor2=-0.0766; - - string description = "If hit wih largest R is less than this, then a "; - description += "fake point will be added when fitting the parameters "; - description += "for the track candidate in the 'FDCCathodes' factory. "; - description += "The point will be on the beamline (x,y) = (0,0) and "; - description += "at a z location determined from the geometry center of "; - description += "target (via DGeometry::GetTargetZ()"; - MAX_R_VERTEX_LIMIT = 50.0; - gPARMS->SetDefaultParameter("TRKFIND:MAX_R_VERTEX_LIMIT", MAX_R_VERTEX_LIMIT, description); - - MATCHING_PHI_CUT=10.0; - gPARMS->SetDefaultParameter("TRKFIND:MATCHING_PHI_CUT", MATCHING_PHI_CUT); + BEAM_VAR=1.; + gPARMS->SetDefaultParameter("TRKFIND:BEAM_VAR",BEAM_VAR); if(DEBUG_HISTS) { dapp->Lock(); @@ -163,7 +132,7 @@ jerror_t DTrackCandidate_factory_FDCCathodes::evnt(JEventLoop *loop, uint64_t ev // abort if there are no segments if (segments.size()==0.) return NOERROR; - std::sort(segments.begin(), segments.end(), DTrackCandidate_segment_cmp); + std::stable_sort(segments.begin(), segments.end(), DTrackCandidate_segment_cmp); // Group segments by package vectorpackages[4]; @@ -191,6 +160,7 @@ jerror_t DTrackCandidate_factory_FDCCathodes::evnt(JEventLoop *loop, uint64_t ev for (unsigned int i=0;iis_quadrupled(triplets.size()); vector >quadruplets; for (unsigned int i=0;i >mytracks; - if (quadruplets.size()==1){ - mytracks.push_back(quadruplets[0]); - } - else if (quadruplets.size()>1){ - // Because segments could have been added to the triplets on either end, - // we need to check for clones - vectoris_clone(quadruplets.size()); - for (unsigned int i=0;iq; + // create a guess for rc and add hits for (unsigned int m=0;mrc; + xc+=mytracks[i][m]->xc; + yc+=mytracks[i][m]->yc; + tanl+=mytracks[i][m]->tanl; for (unsigned int n=0;nhits.size();n++){ const DFDCPseudo *hit=mytracks[i][m]->hits[n]; fit.AddHit(hit); @@ -299,13 +242,12 @@ jerror_t DTrackCandidate_factory_FDCCathodes::evnt(JEventLoop *loop, uint64_t ev } } } - rc/=double(mytracks[i].size()); - // Fake point at origin - bool use_fake_point=false; - if (max_rhits[0]; double Bz=fabs(bfield->GetBz(myhit->xy.X(),myhit->xy.Y(),myhit->wire->origin.z())); double p=0.003*fit.r0*Bz/cos(atan(fit.tanl)); - if (p>10.){//... try alternate circle fit + + // Prune the fake hit at the origin in case we need to use an alternate + // fit + fit.PruneHit(0); + if (p>10.){//... try alternate circle fit fit.FitCircle(); rc=fit.r0; xc=fit.x0; yc=fit.y0; + } + if (rc<0.5*max_r && max_r<10.0){ + // ... we can also have issues near the beam line: + // Try to fix relatively high momentum tracks in the very forward + // direction that look like low momentum tracks due to small pt. + // Assume that the particle came from the center of the target. + fit.FitTrack_FixedZvertex(TARGET_Z); + tanl=fit.tanl; + rc=fit.r0; + xc=fit.x0; + yc=fit.y0; + fit.FindSenseOfRotation(); + q=FactorForSenseOfRotation*fit.h; } - - } - // Try to fix relatively high momentum tracks in the very forward - // direction that look like low momentum tracks due to small pt. - // Assume that the particle came from the center of the target. - - if (rc<0.5*max_r && max_r<10.0){ - if (use_fake_point) fit.PruneHit(0); - fit.FitTrack_FixedZvertex(TARGET_Z); - tanl=fit.tanl; - rc=fit.r0; - xc=fit.x0; - yc=fit.y0; - fit.FindSenseOfRotation(); - q=FactorForSenseOfRotation*fit.h; - } - // Create new track, starting with the most upstream segment DTrackCandidate *track = new DTrackCandidate; //circle fit parameters @@ -355,12 +297,6 @@ jerror_t DTrackCandidate_factory_FDCCathodes::evnt(JEventLoop *loop, uint64_t ev DVector3 mom,pos; GetPositionAndMomentum(mytracks[i],pos,mom); - // Empirical correction to the momentum - if (APPLY_MOMENTUM_CORRECTION){ - double p_mag=mom.Mag(); - mom.SetMag(p_mag*(1.+p_factor1/mom.Theta()+p_factor2)); - } - track->chisq=fit.chisq; track->Ndof=fit.ndof; track->setPID((q > 0.0) ? PiPlus : PiMinus); @@ -375,6 +311,7 @@ jerror_t DTrackCandidate_factory_FDCCathodes::evnt(JEventLoop *loop, uint64_t ev } + // Now try to attach stray segments to existing tracks for (unsigned int i=0;i<4;i++){ for (unsigned int k=0;kFill(my_p,circle_center_diff2_min); + match_center_dist2->Fill(p,circle_center_diff2_min); } return match; } @@ -527,7 +458,6 @@ jerror_t DTrackCandidate_factory_FDCCathodes::GetPositionAndMomentum(const DFDCS double z0=segment->z_vertex; twokappa=FactorForSenseOfRotation*segment->q/segment->rc; - one_over_twokappa=1./twokappa; cotl=1./my_tanl; // Useful intermediate variables @@ -662,7 +592,8 @@ void DTrackCandidate_factory_FDCCathodes::LinkSegments(unsigned int pack1, // Try matching to the next package if (packages[pack2].size()>0 && (match2=GetTrackMatch(segment,packages[pack2],match_id))!=NULL){ - + if (is_paired[pack2][match_id]) continue; + pair mypair; mypair.first=segment; mypair.second=match2; @@ -728,19 +659,27 @@ bool DTrackCandidate_factory_FDCCathodes::LinkStraySegment(const DFDCSegment *se // Add the new segment and sort by z segments.push_back(segment); - sort(segments.begin(),segments.end(),DTrackCandidate_segment_cmp_by_z); + stable_sort(segments.begin(),segments.end(),DTrackCandidate_segment_cmp_by_z); // Create fit object and add hits - DHelicalFit fit; - for (unsigned int m=0;mhits.size();k++){ const DFDCPseudo *hit=segments[m]->hits[k]; fit.AddHit(hit); + + double r=hit->xy.Mod(); + if (r>max_r){ + max_r=r; + } } } // Redo the helical fit with the additional hits - if (fit.FitTrackRiemann(segment->rc)==NOERROR){ + if (fit.FitTrackRiemann(_data[i]->rc)==NOERROR){ rc=fit.r0; tanl=fit.tanl; xc=fit.x0; @@ -751,21 +690,32 @@ bool DTrackCandidate_factory_FDCCathodes::LinkStraySegment(const DFDCSegment *se const DFDCPseudo *myhit=segments[0]->hits[0]; double Bz=fabs(bfield->GetBz(myhit->xy.X(),myhit->xy.Y(),myhit->wire->origin.z())); double p=0.003*fit.r0*Bz/cos(atan(fit.tanl)); + + // Prune the fake hit at the origin in case we need to use an + // alternate fit + fit.PruneHit(0); if (p>10.){//... try alternate circle fit fit.FitCircle(); rc=fit.r0; xc=fit.x0; yc=fit.y0; + } + if (rc<0.5*max_r && max_r<10.0){ + // ... we can also have issues near the beam line: + // Try to fix relatively high momentum tracks in the very forward + // direction that look like low momentum tracks due to small pt. + // Assume that the particle came from the center of the target. + fit.FitTrack_FixedZvertex(TARGET_Z); + tanl=fit.tanl; + rc=fit.r0; + xc=fit.x0; + yc=fit.y0; + fit.FindSenseOfRotation(); + q=FactorForSenseOfRotation*fit.h; } - + // Get position and momentum just upstream of first hit GetPositionAndMomentum(segments,pos,mom); - - // Empirical correction to the momentum - if (APPLY_MOMENTUM_CORRECTION){ - double p_mag=mom.Mag(); - mom.SetMag(p_mag*(1.+p_factor1/mom.Theta()+p_factor2)); - } _data[i]->chisq=fit.chisq; _data[i]->Ndof=fit.ndof; diff --git a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h index 307dc2fd2..3f69fadff 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h +++ b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h @@ -75,27 +75,21 @@ class DTrackCandidate_factory_FDCCathodes:public JFactory{ const DFDCSegment *segment); bool LinkStraySegment(const DFDCSegment *segment); - double MATCHING_PHI_CUT; - bool DEBUG_HISTS,USE_FDC,APPLY_MOMENTUM_CORRECTION; - double p_factor1,p_factor2; + bool DEBUG_HISTS,USE_FDC; + TH2F *match_dist_fdc,*match_center_dist2; - TH1F *Hcircle_fit_prob; + vectorz_wires; - double endplate_z; - double TARGET_Z; - double MAX_R_VERTEX_LIMIT; - double zpack[4]; + double TARGET_Z,BEAM_VAR; + double FactorForSenseOfRotation; // Fit parameters - double xc,yc,rc,z_vertex,q,phi0,tanl; - - // vector of fit results - vectorfit_results; - + double xc,yc,rc,q,tanl; + // Parameters at the end of the segment double xs,ys,zs; - double p,cosphi,sinphi,twokappa,one_over_twokappa,cotl; + double p,cosphi,sinphi,twokappa,cotl; }; inline double DTrackCandidate_factory_FDCCathodes::Match(double p){ diff --git a/src/libraries/TRACKING/DTrackFitter.cc b/src/libraries/TRACKING/DTrackFitter.cc index 537184a33..ce4ed1857 100644 --- a/src/libraries/TRACKING/DTrackFitter.cc +++ b/src/libraries/TRACKING/DTrackFitter.cc @@ -335,31 +335,11 @@ DTrackFitter::FindHitsAndFitTrack(const DKinematicData &starting_params, DTrackHitSelector::fit_type_t input_type = fit_type==kTimeBased ? DTrackHitSelector::kWireBased:DTrackHitSelector::kHelical; hitselector->GetAllHits(input_type, rt, cdctrackhits, fdcpseudos, this,N); - // If the hit selector found no hits at all on the track, the most - // likely explanation is that the charge of the candidate was wrong, - // especially for stiff forward-going tracks. Try rotating phi by - // 180 degrees, switching the charge, and trying the hit selector again. - - if (fdchits.size()+cdchits.size()==0){ - // Make a temporary copy of the current reference trajectory - DReferenceTrajectory temp_rt = *rt; - - // Swim a reference trajectory with revised paramters - mom.SetPhi(mom.Phi()+M_PI); - q*=-1.; - temp_rt.Swim(pos, mom,q); - if(temp_rt.Nswim_steps<1)return fit_status = kFitFailed; - hitselector->GetAllHits(input_type, &temp_rt, cdctrackhits, fdcpseudos, this,N); - - if (fdchits.size()+cdchits.size()!=0){ - if (DEBUG_LEVEL>0) - _DBG_ << "Switching the charge and phi of the track..." <t-b->t)hit->t_u+a->hit->t_v; double tsum_2=b->hit->t_u+b->hit->t_v; + if (fabs(tsum_1-tsum_2)dE>b->dE); + } return (tsum_1tt); @@ -663,8 +666,6 @@ void DTrackFitterKalmanSIMD::ResetKalmanSIMD(void) mT0=0.,mT0MinimumDriftTime=1e6; - mMinDriftTime=1e6; - mMinDriftID=2000; mVarT0=25.; mCDCInternalStepSize=0.5; @@ -714,14 +715,6 @@ DTrackFitter::fit_status_t DTrackFitterKalmanSIMD::FitTrack(void) if (num_good_cdchits>0){ stable_sort(my_cdchits.begin(),my_cdchits.end(),DKalmanSIMDCDCHit_cmp); - // Find earliest time to use for estimate for T0 - for (unsigned int i=0;itdrifttdrift; - mMinDriftID=1000+i; - } - } - // Look for multiple hits on the same wire for (unsigned int i=0;ihit->wire->ring==my_cdchits[i+1]->hit->wire->ring && @@ -741,14 +734,6 @@ DTrackFitter::fit_status_t DTrackFitterKalmanSIMD::FitTrack(void) if (num_good_fdchits>0){ stable_sort(my_fdchits.begin(),my_fdchits.end(),DKalmanSIMDFDCHit_cmp); - // Find earliest time to use for estimate for T0 - for (unsigned int i=0;itt; - } - } - // Look for multiple hits on the same wire for (unsigned int i=0;ihit->wire->layer==my_fdchits[i+1]->hit->wire->layer && @@ -789,40 +774,29 @@ DTrackFitter::fit_status_t DTrackFitterKalmanSIMD::FitTrack(void) fdc_used_in_fit=vector(my_fdchits.size()); } - // start time and variance - mT0=mMinDriftTime; if (fit_type==kTimeBased){ mT0=input_params.t0(); - if (mT0>mMinDriftTime){ - mT0=mMinDriftTime; - mVarT0=7.5; - } - else{ - switch(input_params.t0_detector()){ - case SYS_TOF: - mVarT0=0.01; - break; - case SYS_CDC: - mVarT0=7.5; - break; - case SYS_FDC: - mVarT0=7.5; - break; - case SYS_BCAL: - mVarT0=0.25; - break; - default: - mVarT0=0.09; - break; - } + switch(input_params.t0_detector()){ + case SYS_TOF: + mVarT0=0.01; + break; + case SYS_CDC: + mVarT0=7.5; + break; + case SYS_FDC: + mVarT0=7.5; + break; + case SYS_BCAL: + mVarT0=0.25; + break; + default: + mVarT0=0.09; + break; } - - // _DBG_ << SystemName(input_params.t0_detector()) << " " << mT0 <=0){ - DVector2 origin=my_cdchits[id]->origin; - DVector2 dir=my_cdchits[id]->dir; - DVector2 wire_xy=origin+(forward_traj[m].z-my_cdchits[id]->z0wire)*dir; - DVector2 my_xy(forward_traj[m].S(state_x),forward_traj[m].S(state_y)); - doca2=(wire_xy-my_xy).Mod2(); - - if (doca2>old_doca2){ - if (id==min_id){ - double tcorr=1.18; // not sure why needed.. - mT0MinimumDriftTime=my_cdchits[id]->tdrift-old_time+tcorr; - // _DBG_ << "T0 = " << mT0MinimumDriftTime << endl; - break; - } - doca2=1e6; - id--; - } - } - old_doca2=doca2; - old_time=forward_traj[m].t*TIME_UNIT_CONVERSION; - } - } - if (DEBUG_LEVEL>20) { cout << "--- Forward cdc trajectory ---" <=0){ - origin=my_cdchits[id]->origin; - dir=my_cdchits[id]->dir; - DVector2 wire_xy=origin+(central_traj[m].S(state_z)-my_cdchits[id]->z0wire)*dir; - DVector2 my_xy=central_traj[m].xy; - doca2=(wire_xy-my_xy).Mod2(); - - if (doca2>old_doca2){ - if (id==min_id){ - double tcorr=1.18; // not sure why needed.. - mT0MinimumDriftTime=my_cdchits[id]->tdrift-old_time+tcorr; - //_DBG_ << "T0 = " << mT0MinimumDriftTime << endl; - break; - } - doca2=1e6; - id--; - } - } - old_doca2=doca2; - old_time=central_traj[m].t*TIME_UNIT_CONVERSION; - } - } - - if (DEBUG_LEVEL>20) { cout << "---------" << central_traj.size() <<" entries------" <z)GetField(forward_traj[m].S(state_x),forward_traj[m].S(state_y), z,Bx,By,Bz); @@ -1971,59 +1879,6 @@ jerror_t DTrackFitterKalmanSIMD::SetReferenceTrajectory(DMatrix5x1 &S){ } } - // Find estimate for t0 using smallest drift time - if (fit_type==kWireBased){ - if (mMinDriftID<1000){ - mT0Detector=SYS_FDC; - bool found_minimum=false; - for (unsigned int m=0;m0){ - unsigned int first_hit=forward_traj[m].h_id-1; - for (unsigned int n=0;nt-forward_traj[m].t*TIME_UNIT_CONVERSION+tcorr; - //_DBG_ << "T0 = " << mT0MinimumDriftTime << endl; - found_minimum=true; - break; - } - } - } - } - } - else if (my_cdchits.size()>0 && mMinDriftID>=1000){ - mT0Detector=SYS_CDC; - int id=my_cdchits.size()-1; - double old_time=0.,doca2=0.,old_doca2=1e6; - int min_id=mMinDriftID-1000; - for (unsigned int m=0;m=0){ - DVector2 origin=my_cdchits[id]->origin; - DVector2 dir=my_cdchits[id]->dir; - DVector2 wire_xy=origin+(forward_traj[m].z-my_cdchits[id]->z0wire)*dir; - DVector2 my_xy(forward_traj[m].S(state_x),forward_traj[m].S(state_y)); - doca2=(wire_xy-my_xy).Mod2(); - - if (doca2>old_doca2){ - if (id==min_id){ - double tcorr=1.18; // not sure why needed.. - mT0MinimumDriftTime=my_cdchits[id]->tdrift-old_time+tcorr; - //_DBG_ << "T0 = " << mT0MinimumDriftTime << endl; - break; - } - doca2=1e6; - id--; - } - } - old_doca2=doca2; - old_time=forward_traj[m].t*TIME_UNIT_CONVERSION; - } - } - } - if (DEBUG_LEVEL>20) { cout << "--- Forward fdc trajectory ---" <t-mT0 -forward_traj[k].t*TIME_UNIT_CONVERSION; @@ -5729,7 +5582,7 @@ kalman_error_t DTrackFitterKalmanSIMD::KalmanForwardCDC(double anneal,DMatrix5x1 S+=res*K; } // Mark point on ref trajectory with a hit id for the straw - forward_traj[k].h_id=cdc_index+1; + forward_traj[k].h_id=cdc_index+1000; // Store some updated values related to the hit double scale=(skip_ring)?1.:(1.-H*K); @@ -7102,9 +6955,6 @@ kalman_error_t DTrackFitterKalmanSIMD::ForwardFit(const DMatrix5x1 &S0,const DMa last_forward_pulls.assign(forward_pulls.begin(),forward_pulls.end()); } - // Source for t0 guess - mT0Detector=SYS_CDC; - last_fdc_used_in_fit=fdc_used_in_fit; last_cdc_used_in_fit=cdc_used_in_fit; } //iteration @@ -7117,9 +6967,6 @@ kalman_error_t DTrackFitterKalmanSIMD::ForwardFit(const DMatrix5x1 &S0,const DMa chisq_=chisq_forward; ndf_=last_ndf; - // Source for t0 guess - mT0Detector=SYS_CDC; - // output lists of hits used in the fit and fill pull vector cdchits_used_in_fit.clear(); for (unsigned int m=0;m&cdc_pulls){ DMatrix5x5 Cs=C; DMatrix5x5 A; for (unsigned int m=max-1;m>0;m--){ - if (forward_traj[m].h_id>0){ - unsigned int cdc_index=forward_traj[m].h_id-1; + if (forward_traj[m].h_id>999){ + unsigned int cdc_index=forward_traj[m].h_id-1000; if(cdc_used_in_fit[cdc_index] && my_cdchits[cdc_index]->status == good_hit){ if (DEBUG_LEVEL > 5) { _DBG_ << " Smoothing CDC index " << cdc_index << " ring " << my_cdchits[cdc_index]->hit->wire->ring @@ -9187,6 +9022,8 @@ jerror_t DTrackFitterKalmanSIMD::ExtrapolateForwardToOtherDetectors(){ // Deal with points within fiducial volume of chambers unsigned int fdc_plane=0; + mT0Detector=SYS_NULL; + mT0MinimumDriftTime=1e6; for (int k=intersected_start_counter?index_beyond_start_counter:inner_index;k>=0;k--){ double z=forward_traj[k].z; double t=forward_traj[k].t*TIME_UNIT_CONVERSION; @@ -9197,7 +9034,25 @@ jerror_t DTrackFitterKalmanSIMD::ExtrapolateForwardToOtherDetectors(){ double cosl=cos(atan(tanl)); double pt=cosl/fabs(S(state_q_over_p)); double phi=atan2(S(state_ty),S(state_tx)); - + + // Find estimate for t0 using earliest drift time + if (forward_traj[k].h_id>999){ + unsigned int index=forward_traj[k].h_id-1000; + double dt=my_cdchits[index]->tdrift-t; + if (dt0){ + unsigned int index=forward_traj[k].h_id-1; + double dt=my_fdchits[index]->t-t; + if (dt0){ double ds_theta_ms_sq=3.*fabs(forward_traj[k].Q(state_x,state_x)); @@ -9580,6 +9435,8 @@ jerror_t DTrackFitterKalmanSIMD::ExtrapolateCentralToOtherDetectors(){ // Deal with points within fiducial volume of chambers unsigned int fdc_plane=0; + mT0Detector=SYS_NULL; + mT0MinimumDriftTime=1e6; for (int k=index_beyond_start_counter;k>=0;k--){ S=central_traj[k].S; xy=central_traj[k].xy; @@ -9589,6 +9446,16 @@ jerror_t DTrackFitterKalmanSIMD::ExtrapolateCentralToOtherDetectors(){ double pt=1/fabs(S(state_q_over_pt)); double phi=S(state_phi); + // Find estimate for t0 using earliest drift time + if (central_traj[k].h_id>0){ + unsigned int index=central_traj[k].h_id-1; + double dt=my_cdchits[index]->tdrift-t; + if (dt0){ double ds_theta_ms_sq=3.*fabs(central_traj[k].Q(state_D,state_D)); diff --git a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.h b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.h index b93ced8c1..7fd60f956 100644 --- a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.h +++ b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.h @@ -28,11 +28,6 @@ #define BIG 1.0e8 #define EPS2 1.e-4 #define EPS3 1.e-2 -#define BEAM_RADIUS 0.1 -#define MAX_ITER 25 -#define MAX_CHI2 1e16 -#define CDC_BACKWARD_STEP_SIZE 0.5 -#define NUM_ITER 10 #define Z_MIN -100. #define Z_MAX 370.0 #define R_MAX 65.0 @@ -45,18 +40,14 @@ // divisions by the speed of light #define TIME_UNIT_CONVERSION 3.33564095198152014e-02 #define ONE_OVER_C TIME_UNIT_CONVERSION -#define CDC_DRIFT_SPEED 55e-4 -#define VAR_S 0.09 #define Q_OVER_PT_MAX 100. // 10 MeV/c -#define MAX_PATH_LENGTH 500. #define TAN_MAX 10. - +#define MAX_CHI2 1e6 #define MINIMUM_HIT_FRACTION 0.5 #define DELTA_R 1.0 // distance in r to extend the trajectory beyond the last point -#define CDC_VARIANCE 0.0001 #define FDC_CATHODE_VARIANCE 0.000225 #define FDC_ANODE_VARIANCE 0.000225 @@ -64,10 +55,6 @@ #define ONE_SIXTH 0.16666666666666667 #define TWO_THIRDS 0.66666666666666667 -#define CHISQ_DIFF_CUT 20. -#define MAX_DEDX 40. -#define MIN_ITER 2 -#define MIN_CDC_ITER 0 #define MIN_FDC_HITS 2 #define MIN_CDC_HITS 2 @@ -76,8 +63,6 @@ #define DE_PER_STEP 0.001 // in GeV #define BFIELD_FRAC 0.0001 #define MIN_STEP_SIZE 0.1 // in cm -#define CDC_INTERNAL_STEP_SIZE 0.15 // in cm -#define FDC_INTERNAL_STEP_SIZE 0.5 // in cm #define ELECTRON_MASS 0.000511 // GeV @@ -447,10 +432,6 @@ class DTrackFitterKalmanSIMD: public DTrackFitter{ double m_ratio_sq; // .. and its square double two_m_e; // twice the electron mass double m_e_sq; // square of electron mass - - // minimum drift time - double mMinDriftTime; - unsigned int mMinDriftID; // Lorentz deflection parameters double LORENTZ_NR_PAR1,LORENTZ_NR_PAR2,LORENTZ_NZ_PAR1,LORENTZ_NZ_PAR2; diff --git a/src/libraries/TRACKING/DTrackHitSelectorALT2.cc b/src/libraries/TRACKING/DTrackHitSelectorALT2.cc index 405bba6a7..b4fa324fa 100644 --- a/src/libraries/TRACKING/DTrackHitSelectorALT2.cc +++ b/src/libraries/TRACKING/DTrackHitSelectorALT2.cc @@ -414,12 +414,12 @@ void DTrackHitSelectorALT2::GetCDCHits(fit_type_t fit_type, const DReferenceTraj old_straw=hit->wire->straw; // Find the DOCA to this wire - double s; + double s=0.; double doca = rt->DistToRT(hit->wire, &s); if(!isfinite(doca)) continue; if(!isfinite(s))continue; - if (s<=0.) continue; + if (s<1.) continue; if (doca>MAX_DOCA) continue; const DReferenceTrajectory::swim_step_t *last_step = rt->GetLastSwimStep(); @@ -630,12 +630,12 @@ void DTrackHitSelectorALT2::GetFDCHits(fit_type_t fit_type, const DReferenceTraj const DFDCPseudo *hit = *iter; // Find the DOCA to this wire - double s; + double s=0.; double doca = rt->DistToRT(hit->wire, &s); if(!isfinite(doca)) continue; if(!isfinite(s))continue; - if (s<=0.) continue; + if (s<1.) continue; if (doca>MAX_DOCA)continue; const DReferenceTrajectory::swim_step_t *last_step = rt->GetLastSwimStep(); diff --git a/src/libraries/TRACKING/DTrackWireBased_factory.cc b/src/libraries/TRACKING/DTrackWireBased_factory.cc index 86ba54b86..23fd673ed 100644 --- a/src/libraries/TRACKING/DTrackWireBased_factory.cc +++ b/src/libraries/TRACKING/DTrackWireBased_factory.cc @@ -438,8 +438,10 @@ jerror_t DTrackWireBased_factory::evnt(JEventLoop *loop, uint64_t eventnumber) rtv.push_back(new DReferenceTrajectory(fitter->GetDMagneticFieldMap())); } DReferenceTrajectory *rt = rtv[num_used_rts]; - if(locNumInitialReferenceTrajectories == rtv.size()) //didn't create a new one + if(locNumInitialReferenceTrajectories == rtv.size()){ //didn't create a new one rt->Reset(); + } + rt->SetDGeometry(geom); rt->q = candidate->charge(); @@ -602,8 +604,7 @@ void DTrackWireBased_factory::DoFit(unsigned int c_id, // and position rt->SetMass(mass); //rt->Swim(candidate->position(),candidate->momentum(),candidate->charge()); - rt->Rsqmax_exterior=61.*61.; - rt->FastSwim(candidate->position(),candidate->momentum(),candidate->charge(),2000.0,0.,345.); + rt->FastSwimForHitSelection(candidate->position(),candidate->momentum(),candidate->charge()); status=fitter->FindHitsAndFitTrack(*candidate,rt,loop,mass,candidate->Ndof+3); if (/*false && */status==DTrackFitter::kFitNotDone){ diff --git a/src/plugins/Analysis/imaging/JEventProcessor_imaging.cc b/src/plugins/Analysis/imaging/JEventProcessor_imaging.cc index ce5482e51..507a4e590 100644 --- a/src/plugins/Analysis/imaging/JEventProcessor_imaging.cc +++ b/src/plugins/Analysis/imaging/JEventProcessor_imaging.cc @@ -45,8 +45,8 @@ jerror_t JEventProcessor_imaging::init(void) { gDirectory->mkdir("Vertexes")->cd(); - TwoTrackXYZ= new TH3I("TwoTrackXYZ","z vs y vs x",200,-10,10, - 200,-10,10,140,30,100); + TwoTrackXYZ= new TH3I("TwoTrackXYZ","z vs y vs x",400,-10,10, + 400,-10,10,140,30,100); TwoTrackXYZ->SetXTitle("x (cm)"); TwoTrackXYZ->SetYTitle("y (cm)"); TwoTrackXYZ->SetZTitle("z (cm)"); @@ -70,6 +70,23 @@ jerror_t JEventProcessor_imaging::init(void) DocaPull=new TH1F("DocaPull","#deltad/#sigma(#deltad)",100,0.,5); TwoTrackProb=new TH1F("TwoTrackProb","vertex probability",100,0,1.); + TwoTrackZFit=new TH1F("TwoTrackZFit","z for r<0.5 cm",1000,0,200); + TwoTrackZFit->SetXTitle("z [cm]"); + + TwoTrackPocaCutFit=new TH2F("TwoTrackPocaCutFit","2track POCA,doca cut",4000,0,400,650,0,65); + TwoTrackPocaCutFit->SetXTitle("z (cm)"); + TwoTrackPocaCutFit->SetYTitle("r (cm)"); + + TwoTrackXYFit_at_65cm=new TH2F("TwoTrackXYFit_at_65cm","y vs x near 65 cm",400,-5,5,400,-5,5); + TwoTrackXYFit_at_65cm->SetXTitle("x [cm]"); + TwoTrackXYFit_at_65cm->SetYTitle("y [cm]"); + + TwoTrackXYZFit= new TH3I("TwoTrackXYZFit","z vs y vs x",400,-10,10, + 400,-10,10,140,30,100); + TwoTrackXYZFit->SetXTitle("x (cm)"); + TwoTrackXYZFit->SetYTitle("y (cm)"); + TwoTrackXYZFit->SetZTitle("z (cm)"); + gDirectory->cd("../"); return NOERROR; @@ -95,6 +112,13 @@ jerror_t JEventProcessor_imaging::brun(JEventLoop *eventLoop, int32_t runnumber) FIT_VERTEX=true; gPARMS->SetDefaultParameter("IMAGING:FIT_VERTEX",FIT_VERTEX, "Turn on/off vertex fitting"); + FIT_CL_CUT=0.01; + gPARMS->SetDefaultParameter("IMAGING:FIT_CL_CUT",FIT_CL_CUT, "CL cut for vertex fit"); + TRACK_CL_CUT=1e-4; + gPARMS->SetDefaultParameter("IMAGING:TRACK_CL_CUT",TRACK_CL_CUT, "CL cut for tracks"); + DOCA_CUT=1.0; + gPARMS->SetDefaultParameter("IMAGING:DOCA_CUT",DOCA_CUT, "Maximum doca between tracks"); + return NOERROR; } @@ -128,11 +152,11 @@ jerror_t JEventProcessor_imaging::evnt(JEventLoop *loop, uint64_t eventnumber) for (unsigned int i=0;iGet_BestTrackingFOM()->Get_TrackTimeBased(); - if (TMath::Prob(track1->chisq,track1->Ndof)>0.01){ + if (TMath::Prob(track1->chisq,track1->Ndof)>TRACK_CL_CUT){ for (unsigned int j=i+1;jGet_BestTrackingFOM()->Get_TrackTimeBased(); - if (TMath::Prob(track2->chisq,track2->Ndof)>0.01){ + if (TMath::Prob(track2->chisq,track2->Ndof)>TRACK_CL_CUT){ // Make sure there are enough DReferenceTrajectory objects unsigned int locNumInitialReferenceTrajectories = rtv.size(); while(rtv.size()<=num_used_rts){ @@ -169,35 +193,31 @@ jerror_t JEventProcessor_imaging::evnt(JEventLoop *loop, uint64_t eventnumber) DKinematicData kd1=*track1,kd2=*track2; rt1->IntersectTracks(rt2,&kd1,&kd2,pos,doca,var_doca,vertex_chi2,FIT_VERTEX); // rt1->IntersectTracks(rt2,NULL,NULL,pos,doca,var_doca,vertex_chi2); - if (FIT_VERTEX){ - TwoTrackChi2->Fill(vertex_chi2); - vertex_prob=TMath::Prob(vertex_chi2,1); - TwoTrackProb->Fill(vertex_prob); - } + TwoTrackDoca->Fill(doca); DocaPull->Fill(doca/sqrt(var_doca)); + if (docaFill(pos.z(),pos.Perp()); + TwoTrackXYZ->Fill(pos.x(),pos.y(),pos.z()); + if (pos.z()>64.5 && pos.z()<65.5){ + TwoTrackXY_at_65cm->Fill(pos.x(),pos.y()); + } + if (pos.Perp()<0.5){ + TwoTrackZ->Fill(pos.z()); + } - if (vertex_prob>0.1) - { - - TwoTrackDoca->Fill(doca); - //DVector3 dir1=kd1.momentum(); - //dir1.SetMag(1.); - //DVector3 dir2=kd2.momentum(); - //dir2.SetMag(1.); - //TwoTrackRelCosTheta->Fill(dir1.Dot(dir2)); - if (FIT_VERTEX?doca<0.1:doca<1.) - { - double phi=pos.Phi(); - if (phi<-M_PI) phi+=2.*M_PI; - if (phi>M_PI) phi-=2.*M_PI; - - TwoTrackPocaCut->Fill(pos.z(),pos.Perp()); - TwoTrackXYZ->Fill(pos.x(),pos.y(),pos.z()); - if (pos.z()>64.5 && pos.z()<65.5){ - TwoTrackXY_at_65cm->Fill(pos.x(),pos.y()); - } - if (pos.Perp()<0.5){ - TwoTrackZ->Fill(pos.z()); + if (FIT_VERTEX){ + TwoTrackChi2->Fill(vertex_chi2); + vertex_prob=TMath::Prob(vertex_chi2,1); + TwoTrackProb->Fill(vertex_prob); + if (vertex_prob>FIT_CL_CUT){ + TwoTrackPocaCutFit->Fill(pos.z(),pos.Perp()); + TwoTrackXYZFit->Fill(pos.x(),pos.y(),pos.z()); + if (pos.z()>64.5 && pos.z()<65.5){ + TwoTrackXYFit_at_65cm->Fill(pos.x(),pos.y()); + } + if (pos.Perp()<0.5){ + TwoTrackZFit->Fill(pos.z()); + } } } } diff --git a/src/plugins/Analysis/imaging/JEventProcessor_imaging.h b/src/plugins/Analysis/imaging/JEventProcessor_imaging.h index 67a426192..3bbe3287a 100644 --- a/src/plugins/Analysis/imaging/JEventProcessor_imaging.h +++ b/src/plugins/Analysis/imaging/JEventProcessor_imaging.h @@ -32,7 +32,10 @@ class JEventProcessor_imaging:public jana::JEventProcessor{ TH3I *TwoTrackXYZ; TH1F *TwoTrackDoca,*TwoTrackZ,*TwoTrackRelCosTheta; - TH2F *TwoTrackPocaCut; + TH1F *TwoTrackZFit; + TH2F *TwoTrackXYFit_at_65cm; + TH3I *TwoTrackXYZFit; + TH2F *TwoTrackPocaCut,*TwoTrackPocaCutFit; TH2F *TwoTrackXY_at_65cm; TH1F *TwoTrackChi2,*TwoTrackProb; TH1F *DocaPull; @@ -44,6 +47,7 @@ class JEventProcessor_imaging:public jana::JEventProcessor{ const DMagneticFieldMap *bfield; bool FIT_VERTEX; + double TRACK_CL_CUT,FIT_CL_CUT,DOCA_CUT; }; #endif // _JEventProcessor_imaging_