Skip to content

Commit

Permalink
Merge pull request #386 from JeffersonLab/TrackPositionFix
Browse files Browse the repository at this point in the history 
Track position fix
  • Loading branch information
@aaust
aaust authored May 30, 2020
2 parents a13503b + 29254a2 commit 6555e85
Show file tree
Hide file tree
Showing 3 changed files with 74 additions and 18 deletions.
40 changes: 36 additions & 4 deletions src/libraries/HDDM/DEventSourceREST.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,7 +250,11 @@ jerror_t DEventSourceREST::GetObjects(JEvent &event, JFactory_base *factory)
DGeometry* locGeometry = dapp->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
double locTargetCenterZ = 0.0;
locGeometry->GetTargetZ(locTargetCenterZ);


map<string, double> beam_vals;
if (locEventLoop->GetCalib("PHOTON_BEAM/beam_spot",beam_vals))
throw JException("Could not load CCDB table: PHOTON_BEAM/beam_spot");

vector<double> locBeamPeriodVector;
if(locEventLoop->GetCalib("PHOTON_BEAM/RF/beam_period", locBeamPeriodVector))
throw JException("Could not load CCDB table: PHOTON_BEAM/RF/beam_period");
Expand All @@ -272,6 +276,21 @@ jerror_t DEventSourceREST::GetObjects(JEvent &event, JFactory_base *factory)
dTargetCenterZMap[locRunNumber] = locTargetCenterZ;
dBeamBunchPeriodMap[locRunNumber] = locBeamBunchPeriod;
dDIRCChannelStatusMap[locRunNumber] = locDIRCChannelStatus;

dBeamCenterMap[locRunNumber].Set(beam_vals["x"],
beam_vals["y"]);
dBeamDirMap[locRunNumber].Set(beam_vals["dxdz"],
beam_vals["dydz"]);
dBeamZ0Map[locRunNumber]=beam_vals["z"];

jout << "Run " << locRunNumber << " beam spot:"
<< " x=" << dBeamCenterMap[locRunNumber].X()
<< " y=" << dBeamCenterMap[locRunNumber].Y()
<< " z=" << dBeamZ0Map[locRunNumber]
<< " dx/dz=" << dBeamDirMap[locRunNumber].X()
<< " dy/dz=" << dBeamDirMap[locRunNumber].Y()
<< endl;

}
UnlockRead();

Expand Down Expand Up @@ -352,7 +371,6 @@ jerror_t DEventSourceREST::GetObjects(JEvent &event, JFactory_base *factory)
dynamic_cast<JFactory<DDetectorMatches>*>(factory));
}


return OBJECT_NOT_AVAILABLE;
}

Expand Down Expand Up @@ -1103,6 +1121,18 @@ jerror_t DEventSourceREST::Extract_DTrackTimeBased(hddm_r::HDDM *record,
if (factory==NULL) {
return OBJECT_NOT_AVAILABLE;
}

int locRunNumber = locEventLoop->GetJEvent().GetRunNumber();
DVector2 locBeamCenter,locBeamDir;
double locBeamZ0=0.;
LockRead();
{
locBeamCenter = dBeamCenterMap[locRunNumber];
locBeamDir = dBeamDirMap[locRunNumber];
locBeamZ0 = dBeamZ0Map[locRunNumber];
}
UnlockRead();

string tag = (factory->Tag())? factory->Tag() : "";

vector<DTrackTimeBased*> data;
Expand Down Expand Up @@ -1151,15 +1181,17 @@ jerror_t DEventSourceREST::Extract_DTrackTimeBased(hddm_r::HDDM *record,

// Convert from cartesian coordinates to the 5x1 state vector corresponding to the tracking error matrix.
double vect[5];
DVector2 beam_pos=locBeamCenter+(track_pos.Z()-locBeamZ0)*locBeamDir;
DVector2 diff(track_pos.X()-beam_pos.X(),track_pos.Y()-beam_pos.Y());
vect[2]=tan(M_PI_2 - track_mom.Theta());
vect[1]=track_mom.Phi();
double sinphi=sin(vect[1]);
double cosphi=cos(vect[1]);
vect[0]=tra->charge()/track_mom.Perp();
vect[4]=track_pos.Z();
vect[3]=track_pos.Perp();
vect[3]=diff.Mod();

if ((track_pos.X() > 0 && sinphi>0) || (track_pos.Y() <0 && cosphi>0) || (track_pos.Y() >0 && cosphi<0) || (track_pos.X() <0 && sinphi<0))
if ((diff.X() > 0 && sinphi>0) || (diff.Y() <0 && cosphi>0) || (diff.Y() >0 && cosphi<0) || (diff.X() <0 && sinphi<0))
vect[3] *= -1.;
tra->setTrackingStateVector(vect[0], vect[1], vect[2], vect[3], vect[4]);

Expand Down
3 changes: 3 additions & 0 deletions src/libraries/HDDM/DEventSourceREST.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,9 @@ class DEventSourceREST:public JEventSource
int dDIRCMaxChannels;
enum dirc_status_state {GOOD, BAD, NOISY};
map<unsigned int, vector<vector<int>>> dDIRCChannelStatusMap; //unsigned int is run number

map<unsigned int, DVector2> dBeamCenterMap,dBeamDirMap;
map<unsigned int, double> dBeamZ0Map;

DFCALShower_factory *dFCALShowerFactory;
DBCALShower_factory_IU *dBCALShowerFactory;
Expand Down
49 changes: 35 additions & 14 deletions src/libraries/TRACKING/DTrackFitterKalmanSIMD.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1110,7 +1110,8 @@ inline void DTrackFitterKalmanSIMD::GetMomentum(DVector3 &mom){

// Return the "vertex" position (position at which track crosses beam line)
inline void DTrackFitterKalmanSIMD::GetPosition(DVector3 &pos){
pos.SetXYZ(x_,y_,z_);
DVector2 beam_pos=beam_center+(z_-beam_z0)*beam_dir;
pos.SetXYZ(x_+beam_pos.X(),y_+beam_pos.Y(),z_);
}

// Add GEM points
Expand Down Expand Up @@ -6748,14 +6749,18 @@ kalman_error_t DTrackFitterKalmanSIMD::ForwardFit(const DMatrix5x1 &S0,const DMa
}
// Extrapolate to the point of closest approach to the beam line
z_=last_z;
if (sqrt(Slast(state_x)*Slast(state_x)+Slast(state_y)*Slast(state_y))
>EPS2){
DVector2 beam_pos=beam_center+(z_-beam_z0)*beam_dir;
double dx=Slast(state_x)-beam_pos.X();
double dy=Slast(state_y)-beam_pos.Y();
bool extrapolated=false;
if (sqrt(dx*dx+dy*dy)>EPS2){
DMatrix5x5 Ctemp=Clast;
DMatrix5x1 Stemp=Slast;
double ztemp=z_;
if (ExtrapolateToVertex(Stemp,Ctemp)==NOERROR){
Clast=Ctemp;
Slast=Stemp;
extrapolated=true;
}
else{
//_DBG_ << endl;
Expand All @@ -6775,9 +6780,11 @@ kalman_error_t DTrackFitterKalmanSIMD::ForwardFit(const DMatrix5x1 &S0,const DMa
q_over_pt_=q_over_p_/cosl;
phi_=atan2(ty_,tx_);
if (FORWARD_PARMS_COV==false){
DVector2 beam_pos=beam_center+(z_-beam_z0)*beam_dir;
double dx=x_-beam_pos.X();
double dy=y_-beam_pos.Y();
if (extrapolated){
beam_pos=beam_center+(z_-beam_z0)*beam_dir;
dx=x_-beam_pos.X();
dy=y_-beam_pos.Y();
}
D_=sqrt(dx*dx+dy*dy)+EPS;
x_ = dx; y_ = dy;
double cosphi=cos(phi_);
Expand Down Expand Up @@ -7058,9 +7065,14 @@ kalman_error_t DTrackFitterKalmanSIMD::ForwardCDCFit(const DMatrix5x1 &S0,const

// Extrapolate to the point of closest approach to the beam line
z_=zlast;
if (sqrt(Slast(state_x)*Slast(state_x)+Slast(state_y)*Slast(state_y))
>EPS2)
if (ExtrapolateToVertex(Slast,Clast)!=NOERROR) return EXTRAPOLATION_FAILED;
DVector2 beam_pos=beam_center+(z_-beam_z0)*beam_dir;
double dx=Slast(state_x)-beam_pos.X();
double dy=Slast(state_y)-beam_pos.Y();
bool extrapolated=false;
if (sqrt(dx*dx+dy*dy)>EPS2){
if (ExtrapolateToVertex(Slast,Clast)!=NOERROR) return EXTRAPOLATION_FAILED;
extrapolated=true;
}

// Final momentum, positions and tangents
x_=Slast(state_x), y_=Slast(state_y);
Expand All @@ -7074,9 +7086,11 @@ kalman_error_t DTrackFitterKalmanSIMD::ForwardCDCFit(const DMatrix5x1 &S0,const
q_over_pt_=q_over_p_/cosl;
phi_=atan2(ty_,tx_);
if (FORWARD_PARMS_COV==false){
DVector2 beam_pos=beam_center+(z_-beam_z0)*beam_dir;
double dx=x_-beam_pos.X();
double dy=y_-beam_pos.Y();
if (extrapolated){
beam_pos=beam_center+(z_-beam_z0)*beam_dir;
dx=x_-beam_pos.X();
dy=y_-beam_pos.Y();
}
D_=sqrt(dx*dx+dy*dy)+EPS;
x_ = dx; y_ = dy;
double cosphi=cos(phi_);
Expand Down Expand Up @@ -7290,8 +7304,13 @@ kalman_error_t DTrackFitterKalmanSIMD::CentralFit(const DVector2 &startpos,
// the extrolation vector.
extrapolations[SYS_BCAL].clear();
}
if (last_pos.Mod()>0.001){ // in cm

// Extrapolate to the point of closest approach to the beam line
DVector2 beam_pos=beam_center+(Sclast(state_z)-beam_z0)*beam_dir;
bool extrapolated=false;
if ((last_pos-beam_pos).Mod()>EPS2){ // in cm
if (ExtrapolateToVertex(last_pos,Sclast,Cclast)!=NOERROR) return EXTRAPOLATION_FAILED;
extrapolated=true;
}

// output lists of hits used in the fit
Expand All @@ -7312,7 +7331,9 @@ kalman_error_t DTrackFitterKalmanSIMD::CentralFit(const DVector2 &startpos,
y_=last_pos.Y();
z_=Sclast(state_z);
// Find the (signed) distance of closest approach to the beam line
DVector2 beam_pos=beam_center+(z_-beam_z0)*beam_dir;
if (extrapolated){
beam_pos=beam_center+(z_-beam_z0)*beam_dir;
}
double dx=x_-beam_pos.X();
double dy=y_-beam_pos.Y();
D_=sqrt(dx*dx+dy*dy)+EPS;
Expand Down

0 comments on commit 6555e85

Please sign in to comment.