From d4a1953f70bd179981d23fdd221e443d907acc37 Mon Sep 17 00:00:00 2001
From: Holly Szumila <hszumila@jlab.org>
Date: Sat, 7 Oct 2017 12:52:12 -0400
Subject: [PATCH] Modify drift chamber code for alignment and Correcting hit
 times for SHMS

1) Add X,Y offsets parameter for each plane, instead of X,Y offset per chamber
2) ModifyTHcDriftChamber::CorrectHitTimes to work with SHMS chambers
3) Add GetReadoutSide method to THcDCHit class.
---
 src/THcDC.cxx                | 131 +++++++++++++++++----------
 src/THcDC.h                  |  18 +++-
 src/THcDCHit.h               |  62 ++++++++++++-
 src/THcDriftChamber.cxx      | 168 ++++++++++++++++++++++++++---------
 src/THcDriftChamberPlane.cxx |   2 +
 src/THcDriftChamberPlane.h   |   5 ++
 6 files changed, 291 insertions(+), 95 deletions(-)

diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index 7768a21..a5c8036 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -59,7 +59,11 @@ THcDC::THcDC(
   fNChamber = NULL;
   fWireOrder = NULL;
   fDriftTimeSign = NULL;
+  fReadoutTB = NULL;
+  fReadoutLR = NULL;
 
+  fXPos = NULL;
+  fYPos = NULL;
   fZPos = NULL;
   fAlphaAngle = NULL;
   fBetaAngle = NULL;
@@ -104,7 +108,7 @@ void THcDC::Setup(const char* name, const char* description)
   } else {
     fHMSStyleChambers = 0;
   }
-
+  //cout<<"HMS Style??\t"<<fHMSStyleChambers<<endl;
   string planenamelist;
   DBRequest list[]={
     {"dc_num_planes",&fNPlanes, kInt},
@@ -260,7 +264,11 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
   delete [] fNChamber;  fNChamber = new Int_t [fNPlanes]; // Which chamber is this plane
   delete [] fWireOrder;  fWireOrder = new Int_t [fNPlanes]; // Wire readout order
   delete [] fDriftTimeSign;  fDriftTimeSign = new Int_t [fNPlanes];
+  delete [] fReadoutLR;  fReadoutLR = new Int_t [fNPlanes];
+  delete [] fReadoutTB;  fReadoutTB = new Int_t [fNPlanes];
 
+  delete [] fXPos;  fXPos = new Double_t [fNPlanes];
+  delete [] fYPos;  fYPos = new Double_t [fNPlanes];
   delete [] fZPos;  fZPos = new Double_t [fNPlanes];
   delete [] fAlphaAngle;  fAlphaAngle = new Double_t [fNPlanes];
   delete [] fBetaAngle;  fBetaAngle = new Double_t [fNPlanes];
@@ -270,6 +278,11 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
   delete [] fPlaneTimeZero;  fPlaneTimeZero = new Double_t [fNPlanes];
   delete [] fSigma;  fSigma = new Double_t [fNPlanes];
 
+  Bool_t optional = true;
+  fReadoutLR = new Int_t[fNPlanes];
+  fReadoutTB = new Int_t[fNPlanes];
+
+
   DBRequest list[]={
     {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
     {"dc_wire_velocity",&fWireVelocity,kDouble},
@@ -288,6 +301,8 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
     {"dc_chamber_planes", fNChamber, kInt, (UInt_t)fNPlanes},
     {"dc_wire_counting", fWireOrder, kInt, (UInt_t)fNPlanes},
     {"dc_drifttime_sign", fDriftTimeSign, kInt, (UInt_t)fNPlanes},
+    {"dc_readoutLR", fReadoutLR, kInt, (UInt_t)fNPlanes, optional},
+    {"dc_readoutTB", fReadoutTB, kInt, (UInt_t)fNPlanes, optional},
 
     {"dc_zpos", fZPos, kDouble, (UInt_t)fNPlanes},
     {"dc_alpha_angle", fAlphaAngle, kDouble, (UInt_t)fNPlanes},
@@ -313,7 +328,26 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
     {"debugtrackprint", &fdebugtrackprint , kInt},
     {0}
   };
+   for(Int_t ip=0; ip<fNPlanes;ip++) {
+    fReadoutLR[ip] = 0.0;
+    fReadoutTB[ip] = 0.0;
+   }
+
   gHcParms->LoadParmValues((DBRequest*)&list,fPrefix);
+
+  //Set the default plane x,y positions to those of the chamber
+   for(Int_t ip=0; ip<fNPlanes;ip++) {
+    fXPos[ip] = fXCenter[GetNChamber(ip+1)-1];
+    fYPos[ip] = fYCenter[GetNChamber(ip+1)-1];
+   }
+
+   //Load the x,y positions of the planes if they exist (overwrites defaults)
+   DBRequest listOpt[]={
+     {"dc_xpos", fXPos, kDouble, (UInt_t)fNPlanes, optional},
+     {"dc_ypos", fYPos, kDouble, (UInt_t)fNPlanes, optional},
+     {0}
+   };
+   gHcParms->LoadParmValues((DBRequest*)&listOpt,fPrefix);
   if(fNTracksMaxFP <= 0) fNTracksMaxFP = 10;
   // if(fNTracksMaxFP > HNRACKS_MAX) fNTracksMaxFP = NHTRACKS_MAX;
   cout << "Plane counts:";
@@ -400,7 +434,11 @@ void THcDC::DeleteArrays()
   delete [] fNChamber;   fNChamber = NULL;
   delete [] fWireOrder;   fWireOrder = NULL;
   delete [] fDriftTimeSign;   fDriftTimeSign = NULL;
+  delete [] fReadoutLR;   fReadoutLR = NULL;
+  delete [] fReadoutTB;   fReadoutTB = NULL;
 
+  delete [] fXPos;   fXPos = NULL;
+  delete [] fYPos;   fYPos = NULL;
   delete [] fZPos;   fZPos = NULL;
   delete [] fAlphaAngle;   fAlphaAngle = NULL;
   delete [] fBetaAngle;   fBetaAngle = NULL;
@@ -855,7 +893,6 @@ void THcDC::TrackFit()
   // }
 
   Double_t dummychi2 = 1.0E4;
-
   for(UInt_t itrack=0;itrack<fNDCTracks;itrack++) {
     //    Double_t chi2 = dummychi2;
     //    Int_t htrack_fit_num = itrack;
@@ -867,20 +904,21 @@ void THcDC::TrackFit()
       THcDCHit* hit=theDCTrack->GetHit(ihit);
       planes[ihit]=hit->GetPlaneNum()-1;
       if(fFixLR) {
-	if(fFixPropagationCorrection) {
-	  coords[ihit] = hit->GetPos()
-	    + theDCTrack->GetHitLR(ihit)*theDCTrack->GetHitDist(ihit);
-	} else {
-	  coords[ihit] = hit->GetPos()
-	    + theDCTrack->GetHitLR(ihit)*hit->GetDist();
-	}
+    	  if(fFixPropagationCorrection) {
+    		  coords[ihit] = hit->GetPos()
+	    		+ theDCTrack->GetHitLR(ihit)*theDCTrack->GetHitDist(ihit);
+    	  } else {
+    		  coords[ihit] = hit->GetPos()
+	    		+ theDCTrack->GetHitLR(ihit)*hit->GetDist();
+    	  }
+
       } else {
-	if(fFixPropagationCorrection) {
-	  coords[ihit] = hit->GetPos()
-	    + hit->GetLR()*theDCTrack->GetHitDist(ihit);
-	} else {
-	  coords[ihit] = hit->GetCoord();
-	}
+    	  if(fFixPropagationCorrection) {
+    		  coords[ihit] = hit->GetPos()
+	    		+ hit->GetLR()*theDCTrack->GetHitDist(ihit);
+    	  } else {
+    		  coords[ihit] = hit->GetCoord();
+    	  }
       }
     }
 
@@ -890,26 +928,26 @@ void THcDC::TrackFit()
       TVectorD TT(NUM_FPRAY);
       TMatrixD AA(NUM_FPRAY,NUM_FPRAY);
       for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) {
-	TT[irayp] = 0.0;
-	for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
-	  TT[irayp] += (coords[ihit]*
-			fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]])
-	    /pow(fSigma[planes[ihit]],2);
-	}
+    	  TT[irayp] = 0.0;
+    	  for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
+    		  TT[irayp] += (coords[ihit]*
+    				  fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]])
+					  /pow(fSigma[planes[ihit]],2);
+    	  }
       }
       for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) {
-	for(Int_t jrayp=0;jrayp<NUM_FPRAY;jrayp++) {
-	  AA[irayp][jrayp] = 0.0;
-	  if(jrayp<irayp) { // Symmetric
-	    AA[irayp][jrayp] = AA[jrayp][irayp];
-	  } else {
-	    for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
-	      AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]*
-		fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/
-		pow(fSigma[planes[ihit]],2);
-	    }
-	  }
-	}
+    	  for(Int_t jrayp=0;jrayp<NUM_FPRAY;jrayp++) {
+    		  AA[irayp][jrayp] = 0.0;
+    		  if(jrayp<irayp) { // Symmetric
+    			  AA[irayp][jrayp] = AA[jrayp][irayp];
+    		  } else {
+    			  for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
+    				  AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]*
+    						  fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/
+							  pow(fSigma[planes[ihit]],2);
+    			  }
+    		  }
+    	  }
       }
 
       // Solve 4x4 equations
@@ -926,25 +964,24 @@ void THcDC::TrackFit()
 
       // Make sure fCoords, fResiduals, and fDoubleResiduals are clear
       for(Int_t iplane=0;iplane < fNPlanes; iplane++) {
-	Double_t coord=0.0;
-	for(Int_t ir=0;ir<NUM_FPRAY;ir++) {
-	  coord += fPlaneCoeffs[iplane][raycoeffmap[ir]]*dray[ir];
-	  // cout << "ir = " << ir << ", dray[ir] = " << dray[ir] << endl;
-	}
-	theDCTrack->SetCoord(iplane,coord);
+    	  Double_t coord=0.0;
+    	  for(Int_t ir=0;ir<NUM_FPRAY;ir++) {
+    		  coord += fPlaneCoeffs[iplane][raycoeffmap[ir]]*dray[ir];
+    		  // cout << "ir = " << ir << ", dray[ir] = " << dray[ir] << endl;
+    	  }
+    	  theDCTrack->SetCoord(iplane,coord);
       }
       // Compute Chi2 and residuals
       chi2 = 0.0;
       for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
-	Double_t residual = coords[ihit] - theDCTrack->GetCoord(planes[ihit]);
-	// cout << "ihit = " << ihit << ", planes[ihit] = " << planes[ihit]
-	//      << ", coords[ihit] = " << coords[ihit]
-	//      << ", theDCTrack->GetCoord(planes[ihit]) = "
-	//      << theDCTrack->GetCoord(planes[ihit]) << endl;
-	theDCTrack->SetResidual(planes[ihit], residual);
-	// cout << "Getting residual = " << theDCTrack->GetResidual(planes[ihit]) << endl;
-	chi2 += pow(residual/fSigma[planes[ihit]],2);
+    	  Double_t residual = coords[ihit] - theDCTrack->GetCoord(planes[ihit]);
+    	  theDCTrack->SetResidual(planes[ihit], residual);
+
+  //  	  double track_coord = theDCTrack->GetCoord(planes[ihit]);
+//cout<<planes[ihit]<<"\t"<<track_coord<<"\t"<<coords[ihit]<<"\t"<<residual<<endl;
+    	  chi2 += pow(residual/fSigma[planes[ihit]],2);
       }
+
       theDCTrack->SetVector(dray[0], dray[1], 0.0, dray[2], dray[3]);
     }
     theDCTrack->SetChisq(chi2);
diff --git a/src/THcDC.h b/src/THcDC.h
index 5bea4f0..744be72 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -48,6 +48,8 @@ public:
   Int_t GetTdcWinMin(Int_t plane) const { return fTdcWinMin[plane-1];}
   Int_t GetTdcWinMax(Int_t plane) const { return fTdcWinMax[plane-1];}
 
+  Double_t GetXPos(Int_t plane) const { return fXPos[plane-1];}
+  Double_t GetYPos(Int_t plane) const { return fYPos[plane-1];}
   Double_t GetZPos(Int_t plane) const { return fZPos[plane-1];}
   Double_t GetAlphaAngle(Int_t plane) const { return fAlphaAngle[plane-1];}
   Double_t GetBetaAngle(Int_t plane) const { return fBetaAngle[plane-1];}
@@ -59,6 +61,9 @@ public:
   Double_t GetSpacePointCriterion(Int_t chamber) const { return fSpace_Point_Criterion[chamber-1];}
   Double_t GetCentralTime(Int_t plane) const { return fCentralTime[plane-1];}
   Int_t GetDriftTimeSign(Int_t plane) const { return fDriftTimeSign[plane-1];}
+  Int_t GetReadoutLR(Int_t plane) const { return fReadoutLR[plane-1];}
+  Int_t GetReadoutTB(Int_t plane) const { return fReadoutTB[plane-1];}
+
 
   Double_t GetPlaneTimeZero(Int_t plane) const { return fPlaneTimeZero[plane-1];}
   Double_t GetSigma(Int_t plane) const { return fSigma[plane-1];}
@@ -67,10 +72,12 @@ public:
   Double_t GetNSperChan() const { return fNSperChan;}
 
   Double_t GetCenter(Int_t plane) const {
-    Int_t chamber = GetNChamber(plane)-1;
     return
-      fXCenter[chamber]*sin(fAlphaAngle[plane-1]) +
-      fYCenter[chamber]*cos(fAlphaAngle[plane-1]);
+      //fXCenter[chamber]*sin(fAlphaAngle[plane-1]) +
+      //fYCenter[chamber]*cos(fAlphaAngle[plane-1]);
+
+      fXPos[plane-1]*sin(fAlphaAngle[plane-1]) +
+      fYPos[plane-1]*cos(fAlphaAngle[plane-1]);
   }
   //  friend class THaScCalib;
 
@@ -140,7 +147,12 @@ protected:
   Int_t* fNChamber;
   Int_t* fWireOrder;
   Int_t* fDriftTimeSign;
+  Int_t* fReadoutTB;
+  Int_t* fReadoutLR;
+
 
+  Double_t* fXPos;
+  Double_t* fYPos;
   Double_t* fZPos;
   Double_t* fAlphaAngle;
   Double_t* fBetaAngle;
diff --git a/src/THcDCHit.h b/src/THcDCHit.h
index fd40831..8a2b0e1 100644
--- a/src/THcDCHit.h
+++ b/src/THcDCHit.h
@@ -43,7 +43,6 @@ public:
 
   THcDriftChamberPlane* GetWirePlane() const { return fWirePlane; }
 
-
   void     SetWire(THcDCWire * wire) { fWire = wire; ConvertTimeToDist(); }
   void     SetRawTime(Int_t time)     { fRawTime = time; }
   void     SetTime(Double_t time)     { fTime = time; }
@@ -57,6 +56,66 @@ public:
   Int_t    GetChamberNum() const { return fWirePlane->GetChamberNum(); }
   void     SetCorrectedStatus(Int_t c) { fCorrected = c; }
 
+  //Set this correctly
+  Int_t GetReadoutSide() {
+	  Int_t pln = fWirePlane->GetPlaneNum();
+	  Int_t tb = fWirePlane->GetReadoutTB();
+	  Int_t wn = fWire->GetNum();
+
+	  //check if x board
+	  if ((pln>=3 && pln<=4) || (pln>=9 && pln<=10)){
+		  if (tb>0){
+			  if (wn < 49){
+				  fReadoutSide = 4;
+			  }
+			  else {
+				  fReadoutSide = 2;
+			  }
+		  }
+		  else {
+			  if (wn < 33){
+				  fReadoutSide = 2;
+			  }
+			  else {
+				  fReadoutSide = 4;
+			  }
+		  }
+	  }
+	  //else is u board
+	  else{
+		  if (tb>0){
+		  	if (wn < 41){
+		  		fReadoutSide = 4;
+		  	}
+		  	else if (wn >= 41 && wn <= 63){
+		  		fReadoutSide = 3;
+		  	}
+		  	else if (wn >=64 && wn <=69){
+		  		fReadoutSide = 1;
+		  	}
+		  	else {
+		  		fReadoutSide = 2;
+		  	}
+		  }
+		  else {
+			  if (wn < 39){
+				  fReadoutSide = 2;
+			  }
+			  else if (wn >=39 && wn<=44){
+				  fReadoutSide = 1;
+			  }
+			  else if (wn>=45 && wn<=67){
+				  fReadoutSide = 3;
+			  }
+			  else {
+				  fReadoutSide = 4;
+			  }
+		  }
+	  }
+	  return fReadoutSide;
+  }
+
+
 protected:
   static const Double_t kBig;  //!
 
@@ -70,6 +129,7 @@ protected:
   Double_t    fdDist;    // uncertainty in fDist (for chi2 calc)
   Double_t    ftrDist;   // (Perpendicular) distance from the track
   Int_t       fCorrected; // Has this hit been corrected?
+  Int_t       fReadoutSide; // Side where wire is read out. 1-4 is T/R/B/L from beam view for new chambers.
 
   THcDriftChamber* fChamber; //! Pointer to parent wire plane
 
diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx
index 17e8434..0f27508 100644
--- a/src/THcDriftChamber.cxx
+++ b/src/THcDriftChamber.cxx
@@ -885,6 +885,7 @@ void THcDriftChamber::SelectSpacePoints()
 
 void THcDriftChamber::CorrectHitTimes()
 {
+	//cout<<"next event"<<endl;
   // Use the rough hit positions in the chambers to correct the drift time
   // for hits in the space points.
 
@@ -895,48 +896,123 @@ void THcDriftChamber::CorrectHitTimes()
   // time offset per card will cancel much of the error caused by this.  The
   // alternative is to check by card, rather than by plane and this is harder.
   //if (fhdebugflagpr) cout << "In correcthittimes fNSpacePoints = " << fNSpacePoints << endl;
+
   for(Int_t isp=0;isp<fNSpacePoints;isp++) {
-    THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
-    Double_t x = sp->GetX();
-    Double_t y = sp->GetY();
-    for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) {
-      THcDCHit* hit = sp->GetHit(ihit);
-      THcDriftChamberPlane* plane=hit->GetWirePlane();
-
-      // How do we know this correction only gets applied once?  Is
-      // it determined that a given hit can only belong to one space point?
-      Double_t time_corr = plane->GetReadoutX() ?
-	y*plane->GetReadoutCorr()/fWireVelocity :
-	x*plane->GetReadoutCorr()/fWireVelocity;
-
-      //cout << fFixPropagationCorrection << "  Correcting hit " << hit << " " << plane->GetPlaneNum() << " " << isp << "/" << ihit << "  " << x << "," << y << " time coor = " << time_corr<< endl;
-      // Fortran ENGINE does not do this check, so hits can get "corrected"
-      // multiple times if they belong to multiple space points.
-      // To reproduce the precise ENGINE behavior, remove this if condition.
-      if(fFixPropagationCorrection==0) { // ENGINE behavior
-	hit->SetTime(hit->GetTime() - plane->GetCentralTime()
-		     + plane->GetDriftTimeSign()*time_corr);
-	hit->ConvertTimeToDist();
-	//      hit->SetCorrectedStatus(1);
-      } else {
-	// New behavior: Save corrected distance with the hit in the space point
-	// so that the same hit can have a different correction depending on
-	// which space point it is in.
-	//
-	// This is a hack now because the converttimetodist method is connected to the hit
-	// so I compute the corrected time and distance, and then restore the original
-	// time and distance.  Can probably add a method to hit that does a conversion on a time
-	// but does not modify the hit data.
-	Double_t time=hit->GetTime();
-	Double_t dist=hit->GetDist();
-	hit->SetTime(time - plane->GetCentralTime()
-		     + plane->GetDriftTimeSign()*time_corr);
-	hit->ConvertTimeToDist();
-	sp->SetHitDist(ihit, hit->GetDist());
-	hit->SetTime(time);	// Restore time
-	hit->SetDist(dist);	// Restore distance
-      }
-    }
+	  THcSpacePoint* sp = (THcSpacePoint*)(*fSpacePoints)[isp];
+	  Double_t x = sp->GetX();
+	  Double_t y = sp->GetY();
+	  for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) {
+		  THcDCHit* hit = sp->GetHit(ihit);
+		  THcDriftChamberPlane* plane=hit->GetWirePlane();
+
+		  // How do we know this correction only gets applied once?  Is
+		  // it determined that a given hit can only belong to one space point?
+		  Double_t time_corr = plane->GetReadoutX() ?
+				 y*plane->GetReadoutCorr()/fWireVelocity :
+				 x*plane->GetReadoutCorr()/fWireVelocity;
+
+		  // This applies the wire velocity correction for new SHMS chambers --hszumila, SEP17
+		  if (abs(plane->GetReadoutLR())==1){
+			  Int_t pln = hit->GetPlaneNum();
+			  Int_t readoutSide = hit->GetReadoutSide();
+
+			  THcDC* fParent;
+			  fParent = (THcDC*) GetParent();
+			  Double_t posn = hit->GetPos();
+			  //The following values are determined from param file as permutations on planes 5 and 10
+			  Int_t readhoriz = plane->GetReadoutLR();
+			  Int_t readvert = plane->GetReadoutTB();
+
+			  //+x is up and +y is beam right!
+			  double alpha = fParent->GetAlphaAngle(pln);
+			  double xc = posn*TMath::Sin(alpha);
+			  double yc = posn*TMath::Cos(alpha);
+
+			  Double_t wireDistance = plane->GetReadoutX() ?
+				  (abs(y-yc))*abs(plane->GetReadoutCorr()) :
+				  (abs(x-xc))*abs(plane->GetReadoutCorr());
+
+			  //Readout side is based off wiring diagrams
+			  switch (readoutSide){
+			  case 1: //readout from top of chamber
+				  if (x>xc){wireDistance = -readvert*wireDistance;}
+				  else{wireDistance = readvert*wireDistance;}
+
+				  break;
+			  case 2://readout from right of chamber
+				  if (y>yc){wireDistance = -readhoriz*wireDistance;}
+				  else{wireDistance = readhoriz*wireDistance;}
+
+				  break;
+			  case 3: //readout from bottom of chamber
+				  if (xc>x){wireDistance= -readvert*wireDistance;}
+				  else{wireDistance = readvert*wireDistance;}
+
+				  break;
+			  case 4: //readout from left of chamber
+				  if(yc>y){wireDistance = -readhoriz*wireDistance;}
+				  else{wireDistance = readhoriz*wireDistance;}
+
+				  break;
+			  default:
+				  wireDistance = 0.0;
+			  }
+
+			  Double_t timeCorrection = wireDistance/fWireVelocity;
+
+			  if(fFixPropagationCorrection==0) { // ENGINE behavior
+				  Double_t time=hit->GetTime();
+				  hit->SetTime(time - timeCorrection);
+		          hit->ConvertTimeToDist();
+			  }
+			  else{
+				  Double_t time=hit->GetTime();
+				  Double_t dist=hit->GetDist();
+
+				  hit->SetTime(time - timeCorrection);
+				  //double usingOldTime = time-time_corr;
+				  //hit->SetTime(time- time_corr);
+
+				  hit->ConvertTimeToDist();
+				  sp->SetHitDist(ihit, hit->GetDist());
+
+				  hit->SetTime(time);	// Restore time
+				  hit->SetDist(dist);	// Restore distance
+			  }
+
+		  }
+		  /////////////////////////////////////////////////////////////
+
+		  //     if (fhdebugflagpr) cout << "Correcting hit " << hit << " " << plane->GetPlaneNum() << " " << isp << "/" << ihit << "  " << x << "," << y << endl;
+		  // Fortran ENGINE does not do this check, so hits can get "corrected"
+		  // multiple times if they belong to multiple space points.
+		  // To reproduce the precise ENGINE behavior, remove this if condition.
+		  else{
+			  if(fFixPropagationCorrection==0) { // ENGINE behavior
+				  hit->SetTime(hit->GetTime() - plane->GetCentralTime()
+					  + plane->GetDriftTimeSign()*time_corr);
+				  hit->ConvertTimeToDist();
+			  //      hit->SetCorrectedStatus(1);
+			  } else {
+			  // New behavior: Save corrected distance with the hit in the space point
+			  // so that the same hit can have a different correction depending on
+			  // which space point it is in.
+			  //
+			  // This is a hack now because the converttimetodist method is connected to the hit
+			  // so I compute the corrected time and distance, and then restore the original
+			  // time and distance.  Can probably add a method to hit that does a conversion on a time
+			  // but does not modify the hit data.
+			  Double_t time=hit->GetTime();
+			  Double_t dist=hit->GetDist();
+			  hit->SetTime(time - plane->GetCentralTime()
+					  + plane->GetDriftTimeSign()*time_corr);
+			  hit->ConvertTimeToDist();
+			  sp->SetHitDist(ihit, hit->GetDist());
+			  hit->SetTime(time);	// Restore time
+			  hit->SetDist(dist);	// Restore distance
+			  }
+		  }
+	  }
   }
 }
 UInt_t THcDriftChamber::Count1Bits(UInt_t x)
@@ -994,6 +1070,7 @@ void THcDriftChamber::LeftRight()
     if(fHMSStyleChambers) {
       Int_t smallAngOK = (hasy1>=0) && (hasy2>=0);
       if(fSmallAngleApprox !=0 && smallAngOK) { // to small Angle L/R for Y,Y' planes
+
 	if(sp->GetHit(hasy2)->GetPos() <=
 	   sp->GetHit(hasy1)->GetPos()) {
 	  plusminusknown[hasy1] = -1;
@@ -1031,8 +1108,8 @@ void THcDriftChamber::LeftRight()
 	  }
 	}
 	nplusminus = 1 << (nhits-npaired);
-      }
-    }
+      }//end if fSmallAngleApprox!=0
+    }//end else SOS style
     if(nhits < 2) {
       if (fdebugstubchisq) cout << "THcDriftChamber::LeftRight: numhits-2 < 0" << endl;
     } else if (nhits == 2) {
@@ -1098,6 +1175,7 @@ void THcDriftChamber::LeftRight()
             sp->SetStub(stub);
 	  }
 	}
+	///////////////
       } else if (nplaneshit >= fNPlanes-2 && fHMSStyleChambers) { // Two planes missing
 	Double_t chi2 = FindStub(nhits, sp,
 				     plane_list, bitpat, plusminus, stub);
@@ -1117,7 +1195,8 @@ void THcDriftChamber::LeftRight()
 	  }
           sp->SetStub(stub);
 	}
-      } else {
+      }//////////
+      else {
 	if (fhdebugflagpr) cout << "Insufficient planes hit in THcDriftChamber::LeftRight()" << bitpat <<endl;
       }
     } // End loop of pm combinations
@@ -1227,6 +1306,7 @@ THcDriftChamber::~THcDriftChamber()
   }
 }
 
+
 //_____________________________________________________________________________
 void THcDriftChamber::DeleteArrays()
 {
diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx
index 68dd0dd..863cecd 100644
--- a/src/THcDriftChamberPlane.cxx
+++ b/src/THcDriftChamberPlane.cxx
@@ -129,6 +129,8 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
   fCenter = fParent->GetCenter(fPlaneNum);
   fCentralTime = fParent->GetCentralTime(fPlaneNum);
   fDriftTimeSign = fParent->GetDriftTimeSign(fPlaneNum);
+  fReadoutLR = fParent->GetReadoutLR(fPlaneNum);
+  fReadoutTB = fParent->GetReadoutTB(fPlaneNum);
 
   fNSperChan = fParent->GetNSperChan();
 
diff --git a/src/THcDriftChamberPlane.h b/src/THcDriftChamberPlane.h
index 805ecce..35d3a68 100644
--- a/src/THcDriftChamberPlane.h
+++ b/src/THcDriftChamberPlane.h
@@ -69,6 +69,9 @@ public:
   Double_t*    GetPlaneCoef() { return fPlaneCoef; }
   THcDriftChamberPlane(); // for ROOT I/O
   Double_t     CalcWireFromPos(Double_t pos);
+  Int_t        GetReadoutLR() const { return fReadoutLR;}
+  Int_t        GetReadoutTB() const { return fReadoutTB;}
+
 protected:
 
   TClonesArray* fParentHitList;
@@ -100,6 +103,8 @@ protected:
   Double_t fReadoutCorr;
   Double_t fCentralTime;
   Int_t fDriftTimeSign;
+  Int_t fReadoutLR;
+  Int_t fReadoutTB;
 
   Double_t fCenter;
 
-- 
GitLab