diff --git a/examples/PARAM/hcana.param b/examples/PARAM/hcana.param
index 646ac747fa0834ee7ced8d01863f7f86467ba47c..6d3295a3a54711f41f723fa9be39a9059221562d 100644
--- a/examples/PARAM/hcana.param
+++ b/examples/PARAM/hcana.param
@@ -19,3 +19,14 @@ hhodo_plane_names = "1x 1y 2x 2y"
# The following were defined in REPLAY.PARAM
h_recon_coeff_filename = 'PARAM/hms_recon_coeff.dat' ;hms optics matrix
s_recon_coeff_filename = 'PARAM/sos_recon_coeff.dat' ;sos optics matrix
+
+# The following are set to zero to replicate historical ENGINE behavior
+# For new analyses they should be set to 1. If not defined here,
+# hcana will default 1, the new and correct behaviour.
+
+# If 1, Let a hit have different L/R assignment for different space points
+# instead of L/R assignment from first sp it appears in.
+hdc_fix_lr = 0
+# If 1, don't do the the propagation along the wire each time the hit
+# appears in a space point. (Which means the correction accumulates)
+hdc_fix_propcorr = 0
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index 4dc3af6f6caa97b84a259c905a7e1771afb87f01..9c8a4455653b4b6afe3c899486c6d4b4f8cdcbf8 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -67,6 +67,11 @@ THcDC::THcDC(
fPlaneTimeZero = NULL;
fSigma = NULL;
+ // These should be set to zero (in a parameter file) in order to
+ // replicate historical ENGINE behavior
+ fFixLR = 1;
+ fFixPropagationCorrection = 1;
+
fDCTracks = new TClonesArray( "THcDCTrack", 20 );
}
@@ -301,6 +306,8 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
{"yt_track_criterion", &fYtTrCriterion, kDouble},
{"xpt_track_criterion", &fXptTrCriterion, kDouble},
{"ypt_track_criterion", &fYptTrCriterion, kDouble},
+ {"dc_fix_lr", &fFixLR, kInt},
+ {"dc_fix_propcorr", &fFixPropagationCorrection, kInt},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
@@ -565,7 +572,7 @@ void THcDC::LinkStubs()
}
}
fNDCTracks=0; // Number of Focal Plane tracks found
- fDCTracks->Clear();
+ fDCTracks->Clear("C");
Double_t stubminx = 999999;
Double_t stubminy = 999999;
Double_t stubminxp = 999999;
@@ -575,6 +582,7 @@ void THcDC::LinkStubs()
if (fDebugDC) cout << "Joined space points = " << fNSp-1 << endl;
if(!fSingleStub) {
for(Int_t isp1=0;isp1<fNSp-1;isp1++) { // isp1 is index/id in total list of space points
+ THcSpacePoint* sp1 = fSp[isp1];
if (fDebugDC) cout << "Loop thru joined space points " << isp1+1<< endl;
Int_t sptracks=0;
// Now make sure this sp is not already used in a sp.
@@ -584,7 +592,7 @@ void THcDC::LinkStubs()
THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(itrack));
for(Int_t isp=0;isp<theDCTrack->GetNSpacePoints();isp++) {
// isp is index into list of space points attached to theDCTrack
- if(theDCTrack->GetSpacePointID(isp) == isp1) {
+ if(theDCTrack->GetSpacePoint(isp) == sp1) {
tryflag=0;
}
}
@@ -592,10 +600,11 @@ void THcDC::LinkStubs()
if(tryflag) { // SP not already part of a track
Int_t newtrack=1;
for(Int_t isp2=isp1+1;isp2<fNSp;isp2++) {
+ THcSpacePoint* sp2=fSp[isp2];
if (fDebugDC) cout << "second Loop space points " << isp2<< endl;
- if(fSp[isp1]->fNChamber!=fSp[isp2]->fNChamber) {
- Double_t *spstub1=fSp[isp1]->GetStubP();
- Double_t *spstub2=fSp[isp2]->GetStubP();
+ if(sp1->fNChamber!=sp2->fNChamber) {
+ Double_t *spstub1=sp1->GetStubP();
+ Double_t *spstub2=sp2->GetStubP();
Double_t dposx = spstub1[0] - spstub2[0];
Double_t dposy = spstub1[1] - spstub2[1];
Double_t dposxp = spstub1[2] - spstub2[2];
@@ -626,8 +635,8 @@ void THcDC::LinkStubs()
sptracks=0; // Number of tracks with this seed
stub_tracks[sptracks++] = fNDCTracks;
THcDCTrack *theDCTrack = new( (*fDCTracks)[fNDCTracks++]) THcDCTrack(fNPlanes);
- theDCTrack->AddSpacePoint(isp1);
- theDCTrack->AddSpacePoint(isp2);
+ theDCTrack->AddSpacePoint(sp1);
+ theDCTrack->AddSpacePoint(sp2);
if (fDebugDC) cout << " # sp pts combined = " << theDCTrack->GetNSpacePoints() << endl;
if (fDebugDC) cout << " combine sp = " << isp1 << " and " << isp2 << endl;
// Now save the X, Y and XP for the two stubs
@@ -656,11 +665,11 @@ void THcDC::LinkStubs()
for(Int_t isp=0;isp<theDCTrack->GetNSpacePoints();isp++) {
// isp is index of space points in theDCTrack
if (fDebugDC) cout << "looping of previous track = " << isp+1 << endl;
- if(fSp[isp2]->fNChamber ==
- fSp[theDCTrack->GetSpacePointID(isp)]->fNChamber) {
+ if(sp2->fNChamber ==
+ theDCTrack->GetSpacePoint(isp)->fNChamber) {
spoint=isp;
}
- if(isp2==theDCTrack->GetSpacePointID(isp)) {
+ if(sp2==theDCTrack->GetSpacePoint(isp)) {
duppoint=1;
}
} // End loop over sp in tracks with isp1
@@ -668,7 +677,7 @@ void THcDC::LinkStubs()
// add this space point to current track(2)
if(!duppoint) {
if(spoint<0) {
- theDCTrack->AddSpacePoint(isp2);
+ theDCTrack->AddSpacePoint(sp2);
} else {
// If there is another point in the same chamber
// in this track create a new track with all the
@@ -680,9 +689,9 @@ void THcDC::LinkStubs()
if (fDebugDC) cout << "loop over theDCtrack # of sp tps = " << theDCTrack->GetNSpacePoints() << endl;
for(Int_t isp=0;isp<theDCTrack->GetNSpacePoints();isp++) {
if(isp!=spoint) {
- newDCTrack->AddSpacePoint(theDCTrack->GetSpacePointID(isp));
+ newDCTrack->AddSpacePoint(theDCTrack->GetSpacePoint(isp));
} else {
- newDCTrack->AddSpacePoint(isp2);
+ newDCTrack->AddSpacePoint(sp2);
} // End check for dup on copy
if (fDebugDC) cout << "newDCtrack # of sp tps = " << newDCTrack->GetNSpacePoints() << endl;
} // End copy of track
@@ -706,7 +715,7 @@ void THcDC::LinkStubs()
if(fNDCTracks<MAXTRACKS) {
// Need some constructed at thingy
THcDCTrack *newDCTrack = new( (*fDCTracks)[fNDCTracks++]) THcDCTrack(fNPlanes);
- newDCTrack->AddSpacePoint(isp);
+ newDCTrack->AddSpacePoint(fSp[isp]);
} else {
if (fDebugDC) cout << "EPIC FAIL 3: Too many tracks found in THcDC::LinkStubs" << endl;
fNDCTracks=0;
@@ -715,22 +724,6 @@ void THcDC::LinkStubs()
}
}
}
- // Add the list of hits on the track to the track.
- // Looks like it adds all hits for all space points to every track
- for(Int_t itrack=0;itrack<fNDCTracks;itrack++) {
- THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(itrack));
- theDCTrack->ClearHits();
- if (fDebugDC) cout << " Looping thru track add hits track = " << itrack+1 << " with " << theDCTrack->GetNSpacePoints() << " space points "<< endl;
- // Hit list in the track should have been cleared when created.
- for(Int_t isp=0;isp<theDCTrack->GetNSpacePoints();isp++) {
- Int_t spind=theDCTrack->GetSpacePointID(isp);
- if (fDebugDC) cout << " add hits to " << itrack+1 << " sp pt = " << spind << " hits = " << fSp[spind]->GetNHits() <<endl;
- for(Int_t ihit=0;ihit<fSp[spind]->GetNHits();ihit++) {
- theDCTrack->AddHit(fSp[spind]->GetHit(ihit));
- }
- }
- }
- ///
///
if (fDebugDC) cout << " End Linkstubs Found " << fNDCTracks << " tracks"<<endl;
}
@@ -758,10 +751,35 @@ void THcDC::TrackFit()
}
Double_t dummychi2 = 1.0E4;
+
for(Int_t itrack=0;itrack<fNDCTracks;itrack++) {
// Double_t chi2 = dummychi2;
// Int_t htrack_fit_num = itrack;
THcDCTrack *theDCTrack = static_cast<THcDCTrack*>( fDCTracks->At(itrack));
+
+ Double_t coords[theDCTrack->GetNHits()];
+ Int_t planes[theDCTrack->GetNHits()];
+ for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
+ 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();
+ }
+ } else {
+ if(fFixPropagationCorrection) {
+ coords[ihit] = hit->GetPos()
+ + hit->GetLR()*theDCTrack->GetHitDist(ihit);
+ } else {
+ coords[ihit] = hit->GetCoord();
+ }
+ }
+ }
+
cout << " Looping thru track = " << itrack+1 << " Hits = " << theDCTrack->GetNHits() << endl;
theDCTrack->SetNFree(theDCTrack->GetNHits() - NUM_FPRAY);
Double_t chi2 = dummychi2;
@@ -771,11 +789,9 @@ void THcDC::TrackFit()
for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) {
TT[irayp] = 0.0;
for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
- THcDCHit* hit=theDCTrack->GetHit(ihit);
- Int_t plane=hit->GetPlaneNum()-1;
- TT[irayp] += (hit->GetCoord()*
- fPlaneCoeffs[plane][raycoeffmap[irayp]])
- /pow(fSigma[plane],2);
+ TT[irayp] += (coords[ihit]*
+ fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]])
+ /pow(fSigma[planes[ihit]],2);
}
}
for(Int_t irayp=0;irayp<NUM_FPRAY;irayp++) {
@@ -785,11 +801,9 @@ void THcDC::TrackFit()
AA[irayp][jrayp] = AA[jrayp][irayp];
} else {
for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
- THcDCHit* hit=theDCTrack->GetHit(ihit);
- Int_t plane=hit->GetPlaneNum()-1;
- AA[irayp][jrayp] += fPlaneCoeffs[plane][raycoeffmap[irayp]]*
- fPlaneCoeffs[plane][raycoeffmap[jrayp]]/
- pow(fSigma[plane],2);
+ AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]*
+ fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/
+ pow(fSigma[planes[ihit]],2);
}
}
}
@@ -818,18 +832,16 @@ void THcDC::TrackFit()
// Compute Chi2 and residuals
chi2 = 0.0;
for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
- THcDCHit* hit=theDCTrack->GetHit(ihit);
- Int_t plane=hit->GetPlaneNum()-1;
- Double_t residual = hit->GetCoord() - theDCTrack->GetCoord(plane);
- theDCTrack->SetResidual(plane, residual);
- chi2 += pow(residual/fSigma[plane],2);
+ Double_t residual = coords[ihit] - theDCTrack->GetCoord(planes[ihit]);
+ theDCTrack->SetResidual(planes[ihit], residual);
+ chi2 += pow(residual/fSigma[planes[ihit]],2);
}
if (fDebugDC) {
cout << "Hit HDC_WIRE_COORD Fit postiion Residual " << endl;
for(Int_t ihit=0;ihit < theDCTrack->GetNHits();ihit++) {
THcDCHit* hit=theDCTrack->GetHit(ihit);
Int_t plane=hit->GetPlaneNum()-1;
- cout << ihit+1 << " " << hit->GetCoord() << " " << theDCTrack->GetCoord(plane) << " " << theDCTrack->GetResidual(plane) << endl;
+ cout << ihit+1 << " " << coords[ihit] << " " << theDCTrack->GetCoord(plane) << " " << theDCTrack->GetResidual(plane) << endl;
}
}
theDCTrack->SetVector(dray[0], dray[1], 0.0, dray[2], dray[3]);
@@ -867,7 +879,24 @@ void THcDC::TrackFit()
THcDCHit* hit=theDCTrack2->GetHit(ihit);
Int_t plane=hit->GetPlaneNum()-1;
Double_t pos = DpsiFun(ray1,plane);
- theDCTrack1->SetDoubleResidual(plane,hit->GetCoord() - pos);
+ Double_t coord;
+ if(fFixLR) {
+ if(fFixPropagationCorrection) {
+ coord = hit->GetPos()
+ + theDCTrack2->GetHitLR(ihit)*theDCTrack2->GetHitDist(ihit);
+ } else {
+ coord = hit->GetPos()
+ + theDCTrack2->GetHitLR(ihit)*hit->GetDist();
+ }
+ } else {
+ if(fFixPropagationCorrection) {
+ coord = hit->GetPos()
+ + hit->GetLR()*theDCTrack2->GetHitDist(ihit);
+ } else {
+ coord = hit->GetCoord();
+ }
+ }
+ theDCTrack1->SetDoubleResidual(plane,coord - pos);
// hdc_dbl_res(pln) = hdc_double_residual(1,pln) for hists
}
itrack=0;
@@ -877,7 +906,24 @@ void THcDC::TrackFit()
THcDCHit* hit=theDCTrack1->GetHit(ihit);
Int_t plane=hit->GetPlaneNum()-1;
Double_t pos = DpsiFun(ray1,plane);
- theDCTrack2->SetDoubleResidual(plane,hit->GetCoord() - pos);
+ Double_t coord;
+ if(fFixLR) {
+ if(fFixPropagationCorrection) {
+ coord = hit->GetPos()
+ + theDCTrack1->GetHitLR(ihit)*theDCTrack1->GetHitDist(ihit);
+ } else {
+ coord = hit->GetPos()
+ + theDCTrack1->GetHitLR(ihit)*hit->GetDist();
+ }
+ } else {
+ if(fFixPropagationCorrection) {
+ coord = hit->GetPos()
+ + hit->GetLR()*theDCTrack1->GetHitDist(ihit);
+ } else {
+ coord = hit->GetCoord();
+ }
+ }
+ theDCTrack2->SetDoubleResidual(plane,coord - pos);
// hdc_dbl_res(pln) = hdc_double_residual(1,pln) for hists
}
}
diff --git a/src/THcDC.h b/src/THcDC.h
index 6a46a9cf33128480760cca2c172fa2c8db1df32c..b7f40fd375ae8958ac623a2c3fa9ba8ba31488a6 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -61,6 +61,7 @@ public:
Double_t GetPlaneTimeZero(Int_t plane) const { return fPlaneTimeZero[plane-1];}
Double_t GetSigma(Int_t plane) const { return fSigma[plane-1];}
+ Int_t GetFixPropagationCorrectionFlag() const {return fFixPropagationCorrection;}
Double_t GetNSperChan() const { return fNSperChan;}
@@ -84,6 +85,12 @@ protected:
Int_t fNPlanes; // Total number of DC planes
char** fPlaneNames;
Int_t fNChambers;
+ Int_t fFixLR; // If 1, allow a given hit to have different LR
+ // for different space points
+ Int_t fFixPropagationCorrection; // If 1, don't reapply (and accumulate) the
+ // propagation along the wire correction for
+ // each space point a hit occurs in. Keep a
+ // separate correction for each space point.
// Per-event data
Int_t fNhits;
diff --git a/src/THcDCHit.h b/src/THcDCHit.h
index e32cc6159bd9590170df39d4ebb15a175aebcec3..c0bbc7b325b5033d90f17ff6bf638c0b1a2a7375 100644
--- a/src/THcDCHit.h
+++ b/src/THcDCHit.h
@@ -48,7 +48,8 @@ public:
void SetRawTime(Int_t time) { fRawTime = time; }
void SetTime(Double_t time) { fTime = time; }
void SetDist(Double_t dist) { fDist = dist; }
- void SetLeftRight(Int_t lr) { fCoord = GetPos() + lr*fDist;}
+ void SetLeftRight(Int_t lr) { fCoord = GetPos() + lr*fDist; fLR=lr;}
+ Int_t GetLR() { return fLR; }
void SetdDist(Double_t ddist) { fdDist = ddist; }
void SetFitDist(Double_t dist) { ftrDist = dist; }
Int_t GetPlaneNum() const { return fWirePlane->GetPlaneNum(); }
@@ -64,6 +65,7 @@ protected:
Double_t fTime; // Time corrected for time offset of wire (s)
THcDriftChamberPlane* fWirePlane; //! Pointer to parent wire plane
Double_t fDist; // (Perpendicular) Drift Distance
+ Int_t fLR; // +1/-1 which side of wire
Double_t fCoord; // Actual coordinate of hit
Double_t fdDist; // uncertainty in fDist (for chi2 calc)
Double_t ftrDist; // (Perpendicular) distance from the track
diff --git a/src/THcDCTrack.cxx b/src/THcDCTrack.cxx
index 1423c49980cadea8db12fccbf51b877e91057498..92503b0d1e11b6f1589c03be176445d8faa8a62e 100644
--- a/src/THcDCTrack.cxx
+++ b/src/THcDCTrack.cxx
@@ -7,6 +7,7 @@
#include "THcDCHit.h"
#include "THcDCTrack.h"
+#include "THcSpacePoint.h"
THcDCTrack::THcDCTrack(Int_t nplanes) : fnSP(0), fNHits(0)
{
fHits.clear();
@@ -15,17 +16,26 @@ THcDCTrack::THcDCTrack(Int_t nplanes) : fnSP(0), fNHits(0)
fDoubleResiduals.resize(nplanes);
}
-void THcDCTrack::AddHit(THcDCHit * hit)
+void THcDCTrack::AddHit(THcDCHit * hit, Double_t dist, Int_t lr)
{
// Add a hit to the track
- fHits.push_back(hit);
+ Hit newhit;
+ newhit.dchit = hit;
+ newhit.distCorr = dist;
+ newhit.lr = lr;
+ fHits.push_back(newhit);
fNHits++;
}
-void THcDCTrack::AddSpacePoint( Int_t spid )
+void THcDCTrack::AddSpacePoint( THcSpacePoint* sp )
{
// Add to list of space points in this track
// Need a check for maximum spacepoints of 10
- fspID[fnSP++] = spid;
+ fSp[fnSP++] = sp;
+ // Copy all the hits from the space point into the track
+ // Will need to also copy the corrected distance and lr information
+ for(Int_t ihit=0;ihit<sp->GetNHits();ihit++) {
+ AddHit(sp->GetHit(ihit),sp->GetHitDist(ihit),sp->GetHitLR(ihit));
+ }
}
void THcDCTrack::Clear( const Option_t* )
diff --git a/src/THcDCTrack.h b/src/THcDCTrack.h
index d13c774f5ec755e7d34bfb42f145893f4abad505..4c263a47451d62f53813586c023383c2b3c71ffe 100644
--- a/src/THcDCTrack.h
+++ b/src/THcDCTrack.h
@@ -22,14 +22,16 @@ public:
THcDCTrack(Int_t nplanes);
virtual ~THcDCTrack() {};
- virtual void AddHit(THcDCHit * hit);
- virtual void AddSpacePoint(Int_t spid);
+ virtual void AddSpacePoint(THcSpacePoint* sp);
//Get and Set Functions
- Int_t* GetSpacePoints() {return fspID;}
+ // Int_t* GetSpacePoints() {return fspID;}
Int_t GetNSpacePoints() const { return fnSP;}
- Int_t GetSpacePointID(Int_t i) const {return fspID[i];}
- THcDCHit* GetHit(Int_t i) const {return fHits[i];}
+ // Int_t GetSpacePointID(Int_t i) const {return fspID[i];}
+ THcSpacePoint* GetSpacePoint(Int_t i) const {return fSp[i];}
+ THcDCHit* GetHit(Int_t i) const {return fHits[i].dchit;}
+ Double_t GetHitDist(Int_t ihit) { return fHits[ihit].distCorr; };
+ Int_t GetHitLR(Int_t ihit) { return fHits[ihit].lr; };
Int_t GetNHits() const {return fNHits;}
Int_t GetNFree() const {return fNfree;}
Double_t GetCoord(Int_t ip) const {return fCoords[ip];}
@@ -58,10 +60,20 @@ public:
Int_t fnSP; /* Number of space points in this track */
/* Should we have a list of pointers to space points
instead of their indices? */
- Int_t fspID[10]; /* List of space points in this track */
+ // Int_t fspID[10]; /* List of space points in this track */
+ THcSpacePoint* fSp[10]; /* List of space points in this track */
+
Int_t fNHits;
Int_t fNfree; /* Number of degrees of freedom */
- std::vector<THcDCHit*> fHits; /* List of hits for this track */
+ struct Hit {
+ // This is the same structure as in THcSpacePoint. Can we not
+ // duplicate this?
+ THcDCHit* dchit;
+ Double_t distCorr;
+ Int_t lr;
+ };
+ std::vector<Hit> fHits; /* List of hits for this track */
+ //std::vector<THcDCHit*> fHits; /* List of hits for this track */
std::vector<Double_t> fCoords; /* Coordinate on each plane */
std::vector<Double_t> fResiduals; /* Residual on each plane */
std::vector<Double_t> fDoubleResiduals; /* Residual on each plane for single stub mode */
@@ -69,6 +81,8 @@ public:
Double_t fXp_fp, fYp_fp;
Double_t fChi2_fp;
+ virtual void AddHit(THcDCHit * hit, Double_t dist, Int_t lr);
+
private:
// Hide copy ctor and op=
THcDCTrack( const THcDCTrack& );
diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx
index 33c1b657184733f523def66a4e1ca725b3e11a4d..deef48f4d707af3964d8768df45a20789e56ba71 100644
--- a/src/THcDriftChamber.cxx
+++ b/src/THcDriftChamber.cxx
@@ -146,6 +146,7 @@ Int_t THcDriftChamber::ReadDatabase( const TDatime& date )
fMinHits = fParent->GetMinHits(fChamberNum);
fMaxHits = fParent->GetMaxHits(fChamberNum);
fMinCombos = fParent->GetMinCombos(fChamberNum);
+ fFixPropagationCorrection = fParent->GetFixPropagationCorrectionFlag();
if (fDebugDriftCh) cout << "Chamber " << fChamberNum << " Min/Max: " << fMinHits << "/" << fMaxHits << endl;
fSpacePointCriterion = fParent->GetSpacePointCriterion(fChamberNum);
@@ -919,13 +920,29 @@ void THcDriftChamber::CorrectHitTimes()
// 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()) {
+ if(fFixPropagationCorrection==0) { // ENGINE behavior
hit->SetTime(hit->GetTime() - plane->GetCentralTime()
+ plane->GetDriftTimeSign()*time_corr);
hit->ConvertTimeToDist();
- //if (fDebugDriftCh) cout << ihit+1 << " " << plane->GetPlaneNum() << " " << plane->GetChamberNum() << " " << hit->GetTime() << " " << hit->GetDist() << " " << plane->GetCentralTime() << " " << plane->GetDriftTimeSign() << " " << time_corr << endl;
- hit->SetCorrectedStatus(1);
- //}
+ // 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
+ }
}
}
}
@@ -1022,7 +1039,7 @@ void THcDriftChamber::LeftRight()
}
}
if (nplaneshit >= fNPlanes-1) {
- Double_t chi2 = FindStub(nhits, sp->GetHitVectorP(),
+ Double_t chi2 = FindStub(nhits, sp,
plane_list, bitpat, plusminus, stub);
//if(debugging)
@@ -1059,7 +1076,7 @@ void THcDriftChamber::LeftRight()
}
}
} else if (nplaneshit >= fNPlanes-2) { // Two planes missing
- Double_t chi2 = FindStub(nhits, sp->GetHitVectorP(),
+ Double_t chi2 = FindStub(nhits, sp,
plane_list, bitpat, plusminus, stub);
//if(debugging)
//if (fDebugDriftCh) cout << "pmloop=" << pmloop << " Chi2=" << chi2 << endl;
@@ -1100,7 +1117,10 @@ void THcDriftChamber::LeftRight()
// Calculate final coordinate based on plusminusbest
// Update the hit positions in the space points
for(Int_t ihit=0; ihit<nhits; ihit++) {
+ // Save left/right status with the hit and in the space point
+ // In THcDC will decide which to used based on fix_lr flag
sp->GetHit(ihit)->SetLeftRight(plusminusbest[ihit]);
+ sp->SetHitLR(ihit, plusminusbest[ihit]);
}
// Stubs are calculated in rotated coordinate system
@@ -1129,7 +1149,7 @@ void THcDriftChamber::LeftRight()
}
// if(fAA3Inv.find(bitpat) != fAAInv.end()) { // Valid hit combination
//_____________________________________________________________________________
-Double_t THcDriftChamber::FindStub(Int_t nhits, const std::vector<THcDCHit*>* hits,
+Double_t THcDriftChamber::FindStub(Int_t nhits, THcSpacePoint *sp,
Int_t* plane_list, UInt_t bitpat,
Int_t* plusminus, Double_t* stub)
{
@@ -1141,7 +1161,8 @@ Double_t THcDriftChamber::FindStub(Int_t nhits, const std::vector<THcDCHit*>* hi
// This isn't right. dpos only goes up to 2.
for(Int_t ihit=0;ihit<nhits; ihit++) {
- dpos[ihit] = (*hits)[ihit]->GetPos() + plusminus[ihit]*(*hits)[ihit]->GetDist()
+ dpos[ihit] = sp->GetHit(ihit)->GetPos()
+ + plusminus[ihit]*sp->GetHit(ihit)->GetDist()
- fPsi0[plane_list[ihit]];
for(Int_t index=0;index<3;index++) {
TT[index]+= dpos[ihit]*fStubCoefs[plane_list[ihit]][index]
diff --git a/src/THcDriftChamber.h b/src/THcDriftChamber.h
index c6f57c24fb86e878edb8c24d11b21737833913ef..d71b94874578832b011541d731528b4e0ee615bb 100644
--- a/src/THcDriftChamber.h
+++ b/src/THcDriftChamber.h
@@ -22,6 +22,7 @@
//class THaScCalib;
class TClonesArray;
+class THcSpacePoint;
class THcDriftChamber : public THaSubDetector {
@@ -78,6 +79,7 @@ protected:
Double_t fWireVelocity;
Int_t fSmallAngleApprox;
Double_t fStubMaxXPDiff;
+ Int_t fFixPropagationCorrection;
Double_t fXCenter;
Double_t fYCenter;
@@ -108,7 +110,7 @@ protected:
void ChooseSingleHit(void);
void SelectSpacePoints(void);
UInt_t Count1Bits(UInt_t x);
- Double_t FindStub(Int_t nhits, const std::vector<THcDCHit*>* hits,
+ Double_t FindStub(Int_t nhits, THcSpacePoint *sp,
Int_t* plane_list, UInt_t bitpat,
Int_t* plusminus, Double_t* stub);
diff --git a/src/THcSpacePoint.h b/src/THcSpacePoint.h
index 4122cb7a62e71a30fdc301518d57b77cc99e4ffc..04307e11f64a7155eef06f6f2a64162d3072e0ab 100644
--- a/src/THcSpacePoint.h
+++ b/src/THcSpacePoint.h
@@ -20,24 +20,49 @@ public:
}
virtual ~THcSpacePoint() {}
+ struct Hit {
+ THcDCHit* dchit;
+ Double_t distCorr; // Drift distance corrected by propagation along wire
+ Int_t lr; // Left right flag (+/- 1)
+ // Should we copy into here the information from hit
+ // this is likely to be often used?
+ };
+
void SetXY(Double_t x, Double_t y) {fX = x; fY = y;};
void Clear(Option_t* opt="") {fNHits=0; fNCombos=0; fHits.clear();};
void AddHit(THcDCHit* hit) {
- fHits.push_back(hit);
+ Hit newhit;
+ newhit.dchit = hit;
+ newhit.distCorr = 0.0;
+ newhit.lr = 0;
+ fHits.push_back(newhit);
fNHits++;
}
Int_t GetNHits() {return fNHits;};
void SetNHits(Int_t nhits) {fNHits = nhits;};
Double_t GetX() {return fX;};
Double_t GetY() {return fY;};
- THcDCHit* GetHit(Int_t ihit) {return fHits[ihit];};
- std::vector<THcDCHit*>* GetHitVectorP() {return &fHits;};
- void ReplaceHit(Int_t ihit, THcDCHit *hit) {fHits[ihit] = hit;};
+ THcDCHit* GetHit(Int_t ihit) {return fHits[ihit].dchit;};
+ // std::vector<THcDCHit*>* GetHitVectorP() {return &fHits;};
+ //std::vector<Hit>* GetHitStuffVectorP() {return &fHits;};
+ void ReplaceHit(Int_t ihit, THcDCHit *hit) {
+ fHits[ihit].dchit = hit;
+ fHits[ihit].distCorr = 0.0;
+ fHits[ihit].lr = 0;
+ };
+ void SetHitDist(Int_t ihit, Double_t dist) {
+ fHits[ihit].distCorr = dist;
+ };
+ void SetHitLR(Int_t ihit, Int_t lr) {
+ fHits[ihit].lr = lr;
+ };
void SetStub(Double_t stub[4]) {
for(Int_t i=0;i<4;i++) {
fStub[i] = stub[i];
}
};
+ Double_t GetHitDist(Int_t ihit) { return fHits[ihit].distCorr; };
+ Int_t GetHitLR(Int_t ihit) { return fHits[ihit].lr; };
Double_t* GetStubP() { return fStub; };
void IncCombos() { fNCombos++; };
void SetCombos(Int_t ncombos) { fNCombos=ncombos; };
@@ -53,7 +78,9 @@ public:
Double_t fY;
Int_t fNHits;
Int_t fNCombos;
- std::vector<THcDCHit*> fHits;
+ // This adds more indirection to getting hit information.
+ std::vector<Hit> fHits;
+ //std::vector<THcDCHit*> fHits;
Double_t fStub[4];
// Should we also have a pointer back to the chamber object