From 6a33dcf33f94f291ae17d0c48e7c4b44529b6a29 Mon Sep 17 00:00:00 2001
From: "Stephen A. Wood" <saw@jlab.org>
Date: Thu, 26 Sep 2013 09:54:41 -0400
Subject: [PATCH] Parameters to optionally correct small ENGINE DC hit handling
errors. {prefix}dc_fix_lr Historically, in the ENGINE, if a hit is
used in multiple space points/stubs, the left/right assignment for
that hit, which is later used in track fitting, is that assigned for
the last stub encountered with that hit. Properly, the left right
assignment should be allowed to be different in different space
points. If this parameter is set to zero (e.g. in hcana.param), then
the ENGINE behavior is used. For new analyses, it should be set to
one. (Which is the default) {prefix}dc_fix_propcorr When a
hit gets put into a stub, the distance of the hit from the
discriminator can then be estimated. In the engine, a correction to
the drift time (and thus drift distance) is applied. However, if
that hit ends up in another stub, the correction will get applied again,
resulting in a over correction. Setting this flag to 1 will give that
hit a different corection for each stub that it is in.
These flags will default to the new "correct" way of handling hits if
the above parameters are not set in a parameter file. Currently, both
flags are set to zero in hcana.param to replicate the ENGINE behavior.
To implement these changes, the propagation correction and L/R information
for each hit is saved in space point and track classes. This information
is still saved in the hit class, but only used if in ENGINE compatibility
mode.
The THcDCTrack class now saves a list of space point pointers instead
of space point indices.
The AddSpacePoint method now also copies all the hit information into
the track object so that THcDC doesn't need to explicitely copy all
the hits.
The FindStub method, which fits a stub track to a space point is passed
the space point rather than a list of hits
---
examples/PARAM/hcana.param | 11 +++
src/THcDC.cxx | 142 ++++++++++++++++++++++++-------------
src/THcDC.h | 7 ++
src/THcDCHit.h | 4 +-
src/THcDCTrack.cxx | 18 +++--
src/THcDCTrack.h | 28 ++++++--
src/THcDriftChamber.cxx | 37 +++++++---
src/THcDriftChamber.h | 4 +-
src/THcSpacePoint.h | 37 ++++++++--
9 files changed, 214 insertions(+), 74 deletions(-)
diff --git a/examples/PARAM/hcana.param b/examples/PARAM/hcana.param
index 646ac74..6d3295a 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 4dc3af6..9c8a445 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 6a46a9c..b7f40fd 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 e32cc61..c0bbc7b 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 1423c49..92503b0 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 d13c774..4c263a4 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 33c1b65..deef48f 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 c6f57c2..d71b948 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 4122cb7..04307e1 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
--
GitLab