From 3857b76caa5ee2bcb5cec9bf0a658db5eace5ce4 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Fri, 4 May 2018 16:10:30 -0400 Subject: [PATCH 01/28] New histograms and make some cuts configurable --- .../imaging/JEventProcessor_imaging.cc | 82 ++++++++++++------- .../imaging/JEventProcessor_imaging.h | 6 +- 2 files changed, 56 insertions(+), 32 deletions(-) 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_ From f5e8947c75afcbe9a3925d520eb2131286a52c61 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Tue, 12 Jun 2018 17:35:40 -0400 Subject: [PATCH 02/28] (1) remove unused #defines and other unused variables. (2) Make BEAM_VAR a parameter and remove MAX_R_VERTEX_LIMIT -- always use fake point at origin when doing initial fit, but remove this fake point if alternate refitting is needed. (3) check that segment is not already used in a triplet when making a new triplet. (4) Remove some unneeded code that never gets called any more. --- .../DTrackCandidate_factory_FDCCathodes.cc | 154 ++++++++---------- 1 file changed, 64 insertions(+), 90 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc index 6db0de27e..9e9aa06bd 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,12 +34,6 @@ jerror_t DTrackCandidate_factory_FDCCathodes::brun(JEventLoop* eventLoop, _DBG_<< "FDC geometry not available!" <SetDefaultParameter("TRKFIND:APPLY_MOMENTUM_CORRECTION",APPLY_MOMENTUM_CORRECTION); p_factor1=1.61*M_PI/180.; - p_factor2=-0.0766; + 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); + BEAM_VAR=1.; + gPARMS->SetDefaultParameter("TRKFIND:BEAM_VAR",BEAM_VAR); MATCHING_PHI_CUT=10.0; gPARMS->SetDefaultParameter("TRKFIND:MATCHING_PHI_CUT", MATCHING_PHI_CUT); @@ -191,6 +174,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;irc; for (unsigned int n=0;nhits.size();n++){ @@ -300,12 +252,7 @@ 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 @@ -375,6 +322,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;k0 && (match2=GetTrackMatch(segment,packages[pack2],match_id))!=NULL){ - + if (is_paired[pack2][match_id]) continue; + pair mypair; mypair.first=segment; mypair.second=match2; @@ -731,16 +680,24 @@ bool DTrackCandidate_factory_FDCCathodes::LinkStraySegment(const DFDCSegment *se sort(segments.begin(),segments.end(),DTrackCandidate_segment_cmp_by_z); // Create fit object and add hits - DHelicalFit fit; + DHelicalFit fit; + // Fake point at origin + fit.AddHitXYZ(0.,0.,TARGET_Z,BEAM_VAR,BEAM_VAR,0.); + double max_r=0.; 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,13 +708,30 @@ 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); From bbab621322c4a5e073dddd82a992db4fb04ef482 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Tue, 12 Jun 2018 17:37:07 -0400 Subject: [PATCH 03/28] unused variable cleanup --- src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h index 307dc2fd2..14f5f1bb7 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h +++ b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h @@ -82,9 +82,8 @@ class DTrackCandidate_factory_FDCCathodes:public JFactory{ 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 From cea679179fb06fb8b66906d319cbd98515ab6fd4 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Tue, 12 Jun 2018 17:40:55 -0400 Subject: [PATCH 04/28] Make new method that tries to flip the direction of a CDC candidate that appears to be directed backwards but start counter data appears to contradict this and use the new method at the beginning when projecting CDC tracks to the endplate. --- .../TRACKING/DTrackCandidate_factory.cc | 119 ++++++++++-------- .../TRACKING/DTrackCandidate_factory.h | 3 +- 2 files changed, 72 insertions(+), 50 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index 833cbac94..9b412a612 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory.cc @@ -286,25 +286,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 +323,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){ @@ -440,7 +447,7 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) trackcandidates.push_back(can); } } - + for (unsigned int j=0;jmomentum(); 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); @@ -3174,3 +3153,45 @@ 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; + + 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; +} diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.h b/src/libraries/TRACKING/DTrackCandidate_factory.h index f46cd188b..405c192ae 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.h +++ b/src/libraries/TRACKING/DTrackCandidate_factory.h @@ -129,7 +129,8 @@ class DTrackCandidate_factory:public JFactory{ const DFDCSegment *segment2); bool MatchMethod12(DTrackCandidate *srccan,vector &forward_matches, int &num_fdc_cands_remaining); - + bool TryToFlipDirection(vector&scihits, + DVector3 &mom,DVector3 &pos) const; private: const DMagneticFieldMap *bfield; From 8d574534848ba533157b6746b638068233b3997b Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Thu, 14 Jun 2018 08:35:34 -0400 Subject: [PATCH 05/28] Remove unused vector --- src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h index 14f5f1bb7..9c727c786 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h +++ b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h @@ -88,10 +88,7 @@ class DTrackCandidate_factory_FDCCathodes:public JFactory{ // Fit parameters double xc,yc,rc,z_vertex,q,phi0,tanl; - - // vector of fit results - vectorfit_results; - + // Parameters at the end of the segment double xs,ys,zs; double p,cosphi,sinphi,twokappa,one_over_twokappa,cotl; From ff68e32c3e632a26bb33ee0c2b4b5bbb36fa3122 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Thu, 14 Jun 2018 08:41:34 -0400 Subject: [PATCH 06/28] Make sure that all of the helix variables (xc,yc,rc,tanl) are initialized before calling the Riemann fit -- the fit almost never fails (I only happened to catch it once in one file from run 30300), but when it does fail it can cause unreproducible results in an multi-threaded environment. --- .../TRACKING/DTrackCandidate_factory_FDCCathodes.cc | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc index 9e9aa06bd..e9b984fb5 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc @@ -238,9 +238,14 @@ jerror_t DTrackCandidate_factory_FDCCathodes::evnt(JEventLoop *loop, uint64_t ev // Fake point at origin fit.AddHitXYZ(0.,0.,TARGET_Z,BEAM_VAR,BEAM_VAR,0.); double max_r=0.; - rc=0.; // create a guess for rc and add hits + rc=0.,xc=0.,yc=0.,tanl=0.; //initialize helix variables + q=mytracks[i][0]->q; + // 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); @@ -251,7 +256,11 @@ jerror_t DTrackCandidate_factory_FDCCathodes::evnt(JEventLoop *loop, uint64_t ev } } } - rc/=double(mytracks[i].size()); + double mysize=double(mytracks[i].size()); + rc/=mysize; + xc/=mysize; + yc/=mysize; + tanl/=mysize; // Do the fit if (fit.FitTrackRiemann(rc)==NOERROR){ From 7002ddc6f74c7c286f85abbe49a8b0ca4bea1393 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Thu, 14 Jun 2018 08:43:19 -0400 Subject: [PATCH 07/28] Implement direction-changing trick in another place. --- src/libraries/TRACKING/DTrackCandidate_factory.cc | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index 9b412a612..15c0d340f 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory.cc @@ -434,8 +434,17 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) 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;n Date: Thu, 14 Jun 2018 10:18:20 -0400 Subject: [PATCH 08/28] Remove some variables that are no longer needed. Remove optional momentum correction that was never in practice used. --- .../DTrackCandidate_factory_FDCCathodes.cc | 39 ++----------------- .../DTrackCandidate_factory_FDCCathodes.h | 12 +++--- 2 files changed, 8 insertions(+), 43 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc index e9b984fb5..142c7b74e 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc @@ -35,12 +35,6 @@ jerror_t DTrackCandidate_factory_FDCCathodes::brun(JEventLoop* eventLoop, USE_FDC=false; } - // Get the position of the CDC downstream endplate from DGeometry - double endplate_dz,endplate_rmin,endplate_rmax; - dgeom->GetCDCEndplate(endplate_z,endplate_dz,endplate_rmin,endplate_rmax); - endplate_z+=endplate_dz; - - JCalibration *jcalib = dapp->GetJCalibration(runnumber); map targetparms; if (jcalib->Get("TARGET/target_parms",targetparms)==false){ @@ -53,17 +47,9 @@ 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; - BEAM_VAR=1.; gPARMS->SetDefaultParameter("TRKFIND:BEAM_VAR",BEAM_VAR); - MATCHING_PHI_CUT=10.0; - gPARMS->SetDefaultParameter("TRKFIND:MATCHING_PHI_CUT", MATCHING_PHI_CUT); - if(DEBUG_HISTS) { dapp->Lock(); match_dist_fdc=(TH2F*)gROOT->FindObject("match_dist_fdc"); @@ -311,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); @@ -349,12 +329,6 @@ jerror_t DTrackCandidate_factory_FDCCathodes::evnt(JEventLoop *loop, uint64_t ev // Get the momentum and position at a specific z position DVector3 mom, pos; GetPositionAndMomentum(segment,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)); - } // Create new track, starting with the current segment DTrackCandidate *track = new DTrackCandidate; @@ -387,6 +361,7 @@ double DTrackCandidate_factory_FDCCathodes::DocaSqToHelix(const DFDCPseudo *hit) double sin2ks=sin(twoks); double cos2ks=cos(twoks); double one_minus_cos2ks=1.-cos2ks; + double one_over_twokappa=1./twokappa; double x=xs+(cosphi*sin2ks-sinphi*one_minus_cos2ks)*one_over_twokappa; double y=ys+(sinphi*sin2ks+cosphi*one_minus_cos2ks)*one_over_twokappa; @@ -406,7 +381,6 @@ DFDCSegment *DTrackCandidate_factory_FDCCathodes::GetTrackMatch(DFDCSegment *seg // Get the position and momentum at the exit of the package for the // current segment GetPositionAndMomentum(segment); - double my_p=p; // Match to the next package double doca2_min=1e6,doca2; @@ -463,7 +437,7 @@ DFDCSegment *DTrackCandidate_factory_FDCCathodes::GetTrackMatch(DFDCSegment *seg } } if (DEBUG_HISTS){ - match_center_dist2->Fill(my_p,circle_center_diff2_min); + match_center_dist2->Fill(p,circle_center_diff2_min); } return match; } @@ -484,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 @@ -693,7 +666,7 @@ bool DTrackCandidate_factory_FDCCathodes::LinkStraySegment(const DFDCSegment *se // Fake point at origin fit.AddHitXYZ(0.,0.,TARGET_Z,BEAM_VAR,BEAM_VAR,0.); double max_r=0.; - for (unsigned int m=0;mhits.size();k++){ const DFDCPseudo *hit=segments[m]->hits[k]; fit.AddHit(hit); @@ -743,12 +716,6 @@ bool DTrackCandidate_factory_FDCCathodes::LinkStraySegment(const DFDCSegment *se // 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 9c727c786..3f69fadff 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h +++ b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.h @@ -75,23 +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,BEAM_VAR; double FactorForSenseOfRotation; // Fit parameters - double xc,yc,rc,z_vertex,q,phi0,tanl; + 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){ From fa0abd2e0b9ed9bb654dcd1f8f021b854f6141cc Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Fri, 15 Jun 2018 17:29:40 -0400 Subject: [PATCH 09/28] Remove code that tried to rotate the direction of the track if no hits were found to be on the track from the hit selector code using the assumption that the charge of the candidate was wrong, which can happen for stiff tracks -- this is problematic for busy events... --- src/libraries/TRACKING/DTrackFitter.cc | 26 +++----------------------- 1 file changed, 3 insertions(+), 23 deletions(-) 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..." < Date: Fri, 15 Jun 2018 17:31:42 -0400 Subject: [PATCH 10/28] Add method to try to match stray FDC segments to each other, after all other attempts have failed, by removing the assumption that the tracks originate from the target. --- .../TRACKING/DTrackCandidate_factory.cc | 188 +++++++++++++++++- .../TRACKING/DTrackCandidate_factory.h | 6 +- 2 files changed, 182 insertions(+), 12 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index 15c0d340f..4ab028ced 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 @@ -606,20 +611,27 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) // Get the segment data vectorsegments; 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,segments[0],fdctrackcandidates, - forward_matches)){ + 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; } // Redo circle fit assuming track is almost straight - else if (MatchMethod10(j,srccan,segments[0],fdctrackcandidates, + 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; } - else { + // 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; + } + else{ // Not much we can do here -- add to the final list of candidates DTrackCandidate *can = new DTrackCandidate; // circle parameters @@ -629,7 +641,7 @@ 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++){ @@ -1209,7 +1221,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." <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." <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." <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." <1){ + if (DEBUG_LEVEL>0){ _DBG_ << "Match method #9 succeeded..." << endl; } @@ -2714,7 +2726,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; @@ -3081,6 +3093,160 @@ 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); + } + + fit1.FitCircleRiemann(srccan->rc); + + // 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 + fit1.FitLineRiemann(); + + // 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 forcing circle + // to go through (0,0) + DHelicalFit fit2; + for (unsigned int n=0;nhits.size();n++){ + const DFDCPseudo *hit=segments2[0]->hits[n]; + fit2.AddHit(hit); + } + + fit2.FitCircleRiemann(segments2[0]->rc); + + // 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 + fit2.FitLineRiemann(); + + // 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, diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.h b/src/libraries/TRACKING/DTrackCandidate_factory.h index 405c192ae..9b747d497 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.h +++ b/src/libraries/TRACKING/DTrackCandidate_factory.h @@ -128,7 +128,11 @@ 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; From 3e9d8ad03b02b1fd4cafff6a4a9e8258f66e609f Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Sat, 16 Jun 2018 15:39:44 -0400 Subject: [PATCH 11/28] Add version of FastSwim intended to be used with the hit selector. --- .../TRACKING/DReferenceTrajectory.cc | 116 +++++++++++++++++- src/libraries/TRACKING/DReferenceTrajectory.h | 2 +- 2 files changed, 114 insertions(+), 4 deletions(-) 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, From 155a791286bcfd76fffae21f88ad7c681768bdf4 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Sat, 16 Jun 2018 15:40:54 -0400 Subject: [PATCH 12/28] use new FastSwim routine for hit selector --- src/libraries/TRACKING/DTrackWireBased_factory.cc | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) 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){ From 96f14a8a82ee2e5888165614f84530da65b19770 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Sat, 23 Jun 2018 14:30:34 -0400 Subject: [PATCH 13/28] Add new method for matching stray segments together --- .../TRACKING/DTrackCandidate_factory.cc | 107 ++++++++++++------ .../TRACKING/DTrackCandidate_factory.h | 2 + 2 files changed, 77 insertions(+), 32 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index 4ab028ced..4acb723cf 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory.cc @@ -601,38 +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); - 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; - } - // 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; - } - // 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; + 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; @@ -644,8 +642,8 @@ jerror_t DTrackCandidate_factory::evnt(JEventLoop *loop, uint64_t eventnumber) 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); } @@ -3370,3 +3368,48 @@ bool DTrackCandidate_factory::TryToFlipDirection(vector&schits, 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 9b747d497..972ea4465 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.h +++ b/src/libraries/TRACKING/DTrackCandidate_factory.h @@ -135,6 +135,8 @@ class DTrackCandidate_factory:public JFactory{ 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; From 2d8199a013c3760cb2e5c5d052fe03a9e1118bd0 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Mon, 25 Jun 2018 09:32:52 -0400 Subject: [PATCH 14/28] more error checking --- .../TRACKING/DTrackCandidate_factory.cc | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index 4acb723cf..b6f160151 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory.cc @@ -2233,7 +2233,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; @@ -2432,7 +2432,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; @@ -2469,7 +2469,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; @@ -3113,7 +3113,7 @@ bool DTrackCandidate_factory::MatchMethod13(unsigned int src_index, fit1.AddHit(hit); } - fit1.FitCircleRiemann(srccan->rc); + if (fit1.FitCircleRiemann(srccan->rc)!=NOERROR) return false; // Guess for theta and z from input candidates double theta=srccan->momentum().Theta(); @@ -3121,7 +3121,7 @@ bool DTrackCandidate_factory::MatchMethod13(unsigned int src_index, fit1.z_vertex=srccan->position().Z(); // Redo line fit - fit1.FitLineRiemann(); + 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(); @@ -3144,15 +3144,14 @@ bool DTrackCandidate_factory::MatchMethod13(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]; fit2.AddHit(hit); } - fit2.FitCircleRiemann(segments2[0]->rc); + if (fit2.FitCircleRiemann(segments2[0]->rc)!=NOERROR) continue; // Guess for theta and z from input candidates theta=can2->momentum().Theta(); @@ -3160,7 +3159,7 @@ bool DTrackCandidate_factory::MatchMethod13(unsigned int src_index, fit2.z_vertex=srccan->position().Z(); // Redo line fit - fit2.FitLineRiemann(); + if (fit2.FitLineRiemann()!=NOERROR) continue; // Try to match segments by swimming through the field if (MatchMethod11(q,mypos,mymom,fit2,segment,segments2[0])){ From 0dd4961c6b00c4ef5a0b7adc59b6494269325d31 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Wed, 27 Jun 2018 17:43:06 -0400 Subject: [PATCH 15/28] reduce minimum number of CDC hits for candidate to survive --- src/libraries/TRACKING/DTrackCandidate_factory.cc | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index b6f160151..e99620de3 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory.cc @@ -689,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; @@ -2630,8 +2630,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]; @@ -3335,6 +3334,9 @@ DTrackCandidate_factory::GetPositionAndMomentum(double z,DHelicalFit &fit, 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){ From 742772a402530fae1826dc6e15414c2f65b8e0f7 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Wed, 27 Jun 2018 17:44:33 -0400 Subject: [PATCH 16/28] move algorithm that finds minimum flight-time-corrected drift time to extrapolation code and remove some variables and #defines that are no longer needed. --- .../TRACKING/DTrackFitterKalmanSIMD.cc | 246 ++++-------------- .../TRACKING/DTrackFitterKalmanSIMD.h | 19 -- 2 files changed, 55 insertions(+), 210 deletions(-) diff --git a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc index bd078d958..e05293921 100644 --- a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc +++ b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc @@ -663,8 +663,6 @@ void DTrackFitterKalmanSIMD::ResetKalmanSIMD(void) mT0=0.,mT0MinimumDriftTime=1e6; - mMinDriftTime=1e6; - mMinDriftID=2000; mVarT0=25.; mCDCInternalStepSize=0.5; @@ -714,14 +712,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 +731,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 +771,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 +1876,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 +5579,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 +6952,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 +6964,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 +9019,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 +9031,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 +9432,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 +9443,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..f96f5d82e 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,10 +40,7 @@ // 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. @@ -56,7 +48,6 @@ #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; From 376bce9f4ca333227455f0216ab0909cfd6ef720 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Thu, 28 Jun 2018 09:26:51 -0400 Subject: [PATCH 17/28] restore a #define that was too hastily removed... --- src/libraries/TRACKING/DTrackFitterKalmanSIMD.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.h b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.h index f96f5d82e..7fd60f956 100644 --- a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.h +++ b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.h @@ -43,7 +43,7 @@ #define Q_OVER_PT_MAX 100. // 10 MeV/c #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 From c2060a5581723defa322542457529488dfef95b0 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Thu, 28 Jun 2018 09:42:57 -0400 Subject: [PATCH 18/28] make sure s is initalized before calling DistToRT --- src/libraries/TRACKING/DTrackHitSelectorALT2.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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(); From f37990b7efaaf97177bb5c148d1de3b31c1184ef Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Fri, 29 Jun 2018 13:21:49 -0400 Subject: [PATCH 19/28] add another condition for sorting FDC hits --- src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc index e05293921..61bd2dd10 100644 --- a/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc +++ b/src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc @@ -37,6 +37,9 @@ inline bool static DKalmanSIMDFDCHit_cmp(DKalmanSIMDFDCHit_t *a, DKalmanSIMDFDCH if (fabs(a->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); From 1cf8434b5e469a5321c595c3734f646897a46904 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Mon, 2 Jul 2018 17:42:16 -0400 Subject: [PATCH 20/28] move where ndof is computed, add some error checking code and add routine to compute chisq for circle fit that forces the circle to go through (0,0) --- src/libraries/FDC/DFDCSegment_factory.cc | 59 +++++++++++++++++------- src/libraries/FDC/DFDCSegment_factory.h | 9 ++-- 2 files changed, 48 insertions(+), 20 deletions(-) diff --git a/src/libraries/FDC/DFDCSegment_factory.cc b/src/libraries/FDC/DFDCSegment_factory.cc index 6f3708e3c..3067a0dc6 100644 --- a/src/libraries/FDC/DFDCSegment_factory.cc +++ b/src/libraries/FDC/DFDCSegment_factory.cc @@ -116,7 +116,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 @@ -321,7 +321,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 +438,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 +612,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 +621,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; @@ -732,10 +732,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 +749,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; @@ -803,10 +809,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(segment->hits); - LineFit(segment->hits); + if (CircleFit(segment->hits)==NOERROR){ + if (LineFit(segment->hits)==NOERROR){ + chisq=ComputeCircleChiSq(segment->hits); + } + } } - + FillSegmentData(segment); } } @@ -941,11 +950,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); From edbf7eeb52b4fa82d6bee70d6f8d07a8cd72a057 Mon Sep 17 00:00:00 2001 From: Sean Dobbs Date: Tue, 3 Jul 2018 08:10:57 -0400 Subject: [PATCH 21/28] test --- src/libraries/FDC/DFDCSegment_factory.cc | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/libraries/FDC/DFDCSegment_factory.cc b/src/libraries/FDC/DFDCSegment_factory.cc index 6f3708e3c..ac3e7ef22 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]); @@ -716,7 +725,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 @@ -796,7 +805,7 @@ 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 From 586877bd97455d993f0b097ddd3815def0b292d8 Mon Sep 17 00:00:00 2001 From: Sean Dobbs Date: Tue, 3 Jul 2018 16:41:10 -0400 Subject: [PATCH 22/28] Sort CDC hits --- src/libraries/TRACKING/DTrackCandidate_factory.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index e99620de3..6c94386ca 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory.cc @@ -2100,6 +2100,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 From e7f8d18aa4d5e4c1250b009dd997c65e412175a2 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Thu, 5 Jul 2018 11:31:13 -0400 Subject: [PATCH 23/28] Simplify code that adds stray CDC hits to tracks --- .../TRACKING/DTrackCandidate_factory.cc | 43 +++++++++---------- 1 file changed, 20 insertions(+), 23 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index 6c94386ca..14d2899d7 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory.cc @@ -1915,12 +1915,11 @@ void DTrackCandidate_factory::MatchMethod6(DTrackCandidate *can, DVector3 mom=can->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; } } } @@ -1967,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.; @@ -2010,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 } From 3e3d5687dd8c2fe3faf5708ba9fd5b0003f975b4 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Thu, 5 Jul 2018 11:49:36 -0400 Subject: [PATCH 24/28] switch from sort to stable_sort --- .../TRACKING/DTrackCandidate_factory_CDC.cc | 20 +++++++++---------- .../DTrackCandidate_factory_FDCCathodes.cc | 4 ++-- 2 files changed, 12 insertions(+), 12 deletions(-) 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_FDCCathodes.cc b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc index 142c7b74e..b8c06b428 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory_FDCCathodes.cc @@ -132,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]; @@ -659,7 +659,7 @@ 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; From 6e697e40c7b6c17c7036d9de0c07c858db601f42 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Thu, 5 Jul 2018 17:11:43 -0400 Subject: [PATCH 25/28] convert sort -> stable_sort --- .../TRACKING/DTrackCandidate_factory.cc | 42 +++++++++---------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/src/libraries/TRACKING/DTrackCandidate_factory.cc b/src/libraries/TRACKING/DTrackCandidate_factory.cc index 14d2899d7..42ce6e2ba 100644 --- a/src/libraries/TRACKING/DTrackCandidate_factory.cc +++ b/src/libraries/TRACKING/DTrackCandidate_factory.cc @@ -341,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 @@ -432,7 +432,7 @@ 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; @@ -496,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]); @@ -671,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); } @@ -1200,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]; @@ -1210,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 @@ -1309,7 +1309,7 @@ bool DTrackCandidate_factory::MatchMethod1(const DTrackCandidate *fdccan, // Get the cdc hits associated with this cdc candidate vectorcdchits; 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; @@ -1355,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 @@ -1469,7 +1469,7 @@ bool DTrackCandidate_factory::MatchMethod1(const DTrackCandidate *fdccan, // Get hits already linked to this candidate from associated objects vectorcdchits; 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 @@ -1637,7 +1637,7 @@ bool DTrackCandidate_factory::MatchMethod4(const DTrackCandidate *srccan, srccan->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; @@ -1836,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 @@ -2029,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; @@ -2199,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 @@ -2269,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; @@ -2883,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; @@ -2997,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; From 979238e95ca0e6cde2d590f8892c7d83963403a7 Mon Sep 17 00:00:00 2001 From: Sean Dobbs Date: Fri, 6 Jul 2018 16:13:45 -0400 Subject: [PATCH 26/28] Add another stable_sort conversion --- src/libraries/TRACKING/DTrackCandidate_factory_FDC.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From c9faf6c9bb8ff2106a3a9eb98babddf52cf70fc2 Mon Sep 17 00:00:00 2001 From: Sean Dobbs Date: Fri, 6 Jul 2018 16:14:09 -0400 Subject: [PATCH 27/28] Sort TOF Point results by time... for some reason this affects code reproducibility... --- src/libraries/TOF/DTOFPoint_factory.cc | 9 +++++++++ 1 file changed, 9 insertions(+) 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; } From a385ee34f8f1106e21875be7b1c34a6c1177c0c9 Mon Sep 17 00:00:00 2001 From: Simon Taylor Date: Sat, 7 Jul 2018 15:15:57 -0400 Subject: [PATCH 28/28] add code to cut down on segments that appear to have useable hits for the line fit in one FDC plane --- src/libraries/FDC/DFDCSegment_factory.cc | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/libraries/FDC/DFDCSegment_factory.cc b/src/libraries/FDC/DFDCSegment_factory.cc index 5d6446984..2996b9479 100644 --- a/src/libraries/FDC/DFDCSegment_factory.cc +++ b/src/libraries/FDC/DFDCSegment_factory.cc @@ -170,10 +170,12 @@ jerror_t DFDCSegment_factory::RiemannLineFit(vector&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&points, //oldx=XYZ(k,0); //oldy=XYZ(k,1); oldxy=XYZ[k].xy; - } + } Delta=sumv*sumxx-sumx*sumx; double denom=sumv*sumxy-sumy*sumx; @@ -224,6 +226,10 @@ jerror_t DFDCSegment_factory::RiemannLineFit(vector&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; }