diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index 07eb1c6c3534741db4816b3df667b231e25018f9..26163a942b568ceab8363c9f9e9a61b07469e547 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -404,6 +404,7 @@ Int_t THcDC::CoarseTrack( TClonesArray& /* tracks */ )
for(Int_t i=0;i<fNChambers;i++) {
fChambers[i]->FindSpacePoints();
+ fChambers[i]->CorrectHitTimes();
}
ApplyCorrections();
diff --git a/src/THcDC.h b/src/THcDC.h
index d73491b49572d829720744d6f038eda2db06e7f0..aac69980c7ea9708e8b71e14d59ae4aa339acd5f 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -52,6 +52,8 @@ public:
Int_t GetMaxHits(Int_t chamber) const { return fMaxHits[chamber-1];}
Int_t GetMinCombos(Int_t chamber) const { return fMinCombos[chamber-1];}
Double_t GetSpacePointCriterion(Int_t chamber) const { return TMath::Sqrt(fSpace_Point_Criterion2[chamber-1]);}
+ Double_t GetCentralTime(Int_t plane) const { return fCentralTime[plane-1];}
+ Int_t GetDriftTimeSign(Int_t plane) const { return fDriftTimeSign[plane-1];}
Double_t GetPlaneTimeZero(Int_t plane) const { return fPlaneTimeZero[plane-1];}
diff --git a/src/THcDCHit.h b/src/THcDCHit.h
index fb886f6b1138ecea0d8ef53b757a4a15786c7d21..c44e53206e89c85f0332155bd441d6f08447e22f 100644
--- a/src/THcDCHit.h
+++ b/src/THcDCHit.h
@@ -21,6 +21,7 @@ public:
fWire(wire), fRawTime(rawtime), fTime(time), fWirePlane(wp),
fDist(0.0), ftrDist(kBig) {
ConvertTimeToDist();
+ fCorrected = 0;
}
virtual ~THcDCHit() {}
@@ -37,6 +38,7 @@ public:
Double_t GetDist() const { return fDist; }
Double_t GetPos() const { return fWire->GetPos(); } //Position of hit wire
Double_t GetdDist() const { return fdDist; }
+ Int_t GetCorrectedStatus() const { return fCorrected; }
THcDriftChamberPlane* GetWirePlane() const { return fWirePlane; }
@@ -50,6 +52,7 @@ public:
Int_t GetPlaneNum() const { return fWirePlane->GetPlaneNum(); }
Int_t GetPlaneIndex() const { return fWirePlane->GetPlaneIndex(); }
Int_t GetChamberNum() const { return fWirePlane->GetChamberNum(); }
+ void SetCorrectedStatus(Int_t c) { fCorrected = c; }
protected:
static const Double_t kBig; //!
@@ -61,6 +64,7 @@ protected:
Double_t fDist; // (Perpendicular) Drift Distance
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?
THcDriftChamber* fChamber; //! Pointer to parent wire plane
diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx
index 1a77d4d7f230f59aaed0cc5349e44412130cb4dc..7e66401b66981959460d8170c3b1c0141079f508 100644
--- a/src/THcDriftChamber.cxx
+++ b/src/THcDriftChamber.cxx
@@ -119,6 +119,7 @@ Int_t THcDriftChamber::ReadDatabase( const TDatime& date )
prefix[1]='\0';
DBRequest list[]={
{"_remove_sppt_if_one_y_plane",&fRemove_Sppt_If_One_YPlane, kInt},
+ {"dc_wire_velocity", &fWireVelocity, kDouble},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
@@ -240,7 +241,7 @@ Int_t THcDriftChamber::FindSpacePoints( void )
SelectSpacePoints();
if(fNSpacePoints == 0) cout << "SelectSpacePoints() killed SP" << endl;
}
- //cout << fNSpacePoints << " Space Points remain" << endl;
+ cout << fNSpacePoints << " Space Points remain" << endl;
// Add these space points to the total list of space points for the
// the DC package. Do this in THcDC.cxx.
#if 0
@@ -358,12 +359,14 @@ Int_t THcDriftChamber::FindHardSpacePoints()
Combo combos[10*MAX_NUMBER_PAIRS];
for(Int_t ipair1=0;ipair1<ntest_points-1;ipair1++) {
for(Int_t ipair2=ipair1+1;ipair2<ntest_points;ipair2++) {
- Double_t dist2 = pow(pairs[ipair1].x - pairs[ipair2].x,2)
- + pow(pairs[ipair1].y - pairs[ipair2].y,2);
- if(dist2 <= fSpacePointCriterion2) {
- combos[ncombos].pair1 = &pairs[ipair1];
- combos[ncombos].pair2 = &pairs[ipair2];
- ncombos++;
+ if(ncombos < 10*MAX_NUMBER_PAIRS) {
+ Double_t dist2 = pow(pairs[ipair1].x - pairs[ipair2].x,2)
+ + pow(pairs[ipair1].y - pairs[ipair2].y,2);
+ if(dist2 <= fSpacePointCriterion2) {
+ combos[ncombos].pair1 = &pairs[ipair1];
+ combos[ncombos].pair2 = &pairs[ipair2];
+ ncombos++;
+ }
}
}
}
@@ -375,8 +378,8 @@ Int_t THcDriftChamber::FindHardSpacePoints()
hits[2]=combos[icombo].pair2->hit1;
hits[3]=combos[icombo].pair2->hit2;
// Get Average Space point xt, yt
- Double_t xt = combos[icombo].pair1->x + combos[icombo].pair2->x;
- Double_t yt = combos[icombo].pair1->y + combos[icombo].pair2->y;
+ Double_t xt = (combos[icombo].pair1->x + combos[icombo].pair2->x)/2.0;
+ Double_t yt = (combos[icombo].pair1->y + combos[icombo].pair2->y)/2.0;
// Loop over space points
if(fNSpacePoints > 0) {
Int_t add_flag=1;
@@ -418,8 +421,10 @@ Int_t THcDriftChamber::FindHardSpacePoints()
if(iflag[icm]==0) {
fSpacePoints[ispace].hits[fSpacePoints[ispace].nhits++] = hits[icm];
}
- fSpacePoints[ispace].ncombos++;
}
+ fSpacePoints[ispace].ncombos++;
+ // Terminate loop since this combo can only belong to one space point
+ break;
}
}
}// End of loop over existing space points
@@ -709,20 +714,43 @@ void THcDriftChamber::SelectSpacePoints()
fNSpacePoints = sp_count;
}
-/*
-*
-* Now we know rough hit positions in the chambers so we can make
-* wire velocity drift time corrections for each hit in the space point
-*
-* Assume all wires for a plane are read out on the same side (l/r or t/b).
-* If the wire is closer to horizontal, read out left/right. If nearer
-* vertical, assume top/bottom. (Note, this is not always true for the
-* SOS u and v planes. They have 1 card each on the side, but the overall
-* 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.
-*/
-
-
+void THcDriftChamber::CorrectHitTimes()
+{
+ // Use the rough hit positions in the chambers to correct the drift time
+ // for hits in the space points.
+
+ // Assume all wires for a plane are read out on the same side (l/r or t/b).
+ // If the wire is closer to horizontal, read out left/right. If nearer
+ // vertical, assume top/bottom. (Note, this is not always true for the
+ // SOS u and v planes. They have 1 card each on the side, but the overall
+ // 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.
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ Double_t x = fSpacePoints[isp].x;
+ Double_t y = fSpacePoints[isp].y;
+ for(Int_t ihit=0;ihit<fSpacePoints[isp].nhits;ihit++) {
+ THcDCHit* hit = fSpacePoints[isp].hits[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 << "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.
+ if(! hit->GetCorrectedStatus()) {
+ hit->SetTime(hit->GetTime() - plane->GetCentralTime()
+ + plane->GetDriftTimeSign()*time_corr);
+ hit->ConvertTimeToDist();
+ hit->SetCorrectedStatus(1);
+ }
+ }
+ }
+}
//_____________________________________________________________________________
THcDriftChamber::~THcDriftChamber()
{
diff --git a/src/THcDriftChamber.h b/src/THcDriftChamber.h
index 96c293c0fbe4c6a7600a712c3a907a804cc0413f..0e7b6442d9a073aa44db9850c207f1a070de0651 100644
--- a/src/THcDriftChamber.h
+++ b/src/THcDriftChamber.h
@@ -33,6 +33,7 @@ public:
virtual Int_t ApplyCorrections( void );
virtual void ProcessHits( void );
virtual Int_t FindSpacePoints( void ) ;
+ virtual void CorrectHitTimes( void ) ;
virtual void Clear( Option_t* opt="" );
@@ -66,6 +67,7 @@ protected:
Int_t fMaxHits; // Maximum required to do something
Int_t fMinCombos; // Minimum # pairs in a space point
Int_t fRemove_Sppt_If_One_YPlane;
+ Double_t fWireVelocity;
Double_t fXCenter;
Double_t fYCenter;
diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx
index b23cfe21c160216f272ab7c3abd6cf807802dc2d..8c63c5d9b3a03d385d6bc14c9938776c23e385d0 100644
--- a/src/THcDriftChamberPlane.cxx
+++ b/src/THcDriftChamberPlane.cxx
@@ -110,6 +110,8 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
fTdcWinMax = fParent->GetTdcWinMax(fPlaneNum);
fPlaneTimeZero = fParent->GetPlaneTimeZero(fPlaneNum);
fCenter = fParent->GetCenter(fPlaneNum);
+ fCentralTime = fParent->GetCentralTime(fPlaneNum);
+ fDriftTimeSign = fParent->GetDriftTimeSign(fPlaneNum);
fNSperChan = fParent->GetNSperChan();
@@ -133,6 +135,14 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
Double_t hychi = sinalpha*cosgamma;
Double_t hypsi = cosalpha*cosgamma;
+ if(cosalpha <= 0.707) { // x-like wire, need dist from x=0 line
+ fReadoutX = 1;
+ fReadoutCorr = 1/sinalpha;
+ } else {
+ fReadoutX = 0;
+ fReadoutCorr = 1/cosalpha;
+ }
+
Double_t sumsqupsi = hzpsi*hzpsi+hxpsi*hxpsi+hypsi*hypsi;
Double_t sumsquchi = hzchi*hzchi+hxchi*hxchi+hychi*hychi;
Double_t sumcross = hzpsi*hzchi + hxpsi*hxchi + hypsi*hychi;
diff --git a/src/THcDriftChamberPlane.h b/src/THcDriftChamberPlane.h
index 78faf37fefb33635e0c8ac1546c2ea9caf34abe2..e81c9ee5730e37a24169974bbea53fd405c0e0fc 100644
--- a/src/THcDriftChamberPlane.h
+++ b/src/THcDriftChamberPlane.h
@@ -57,6 +57,10 @@ class THcDriftChamberPlane : public THaSubDetector {
Int_t GetPlaneIndex() { return fPlaneIndex; }
Double_t GetXsp() const { return fXsp; }
Double_t GetYsp() const { return fYsp; }
+ Int_t GetReadoutX() { return fReadoutX; }
+ Double_t GetReadoutCorr() { return fReadoutCorr; }
+ Double_t GetCentralTime() { return fCentralTime; }
+ Int_t GetDriftTimeSign() { return fDriftTimeSign; }
protected:
@@ -78,6 +82,11 @@ class THcDriftChamberPlane : public THaSubDetector {
Double_t fXsp;
Double_t fYsp;
+ Int_t fReadoutX;
+ Double_t fReadoutCorr;
+ Double_t fCentralTime;
+ Int_t fDriftTimeSign;
+
Double_t fCenter;
Double_t fNSperChan; /* TDC bin size */