From 7e790fc024febfaaa26853d6831f31757768cadc Mon Sep 17 00:00:00 2001
From: hallc-online <hallconline@gmail.com>
Date: Sat, 17 Mar 2018 07:41:56 -0400
Subject: [PATCH] Modify THcDC and THcDCTrack
Modify THcDCTrack::TrackFit to calculate residuals for each plane
without the plane involved in the track.
Add tree variable residualExclPlane
---
src/THcDC.cxx | 52 ++++++++++++++++++++++++++++++++++++++++++++++
src/THcDC.h | 1 +
src/THcDCTrack.cxx | 1 +
src/THcDCTrack.h | 3 +++
4 files changed, 57 insertions(+)
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index af7e93a..675071a 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -241,6 +241,7 @@ THaAnalysisObject::EStatus THcDC::Init( const TDatime& date )
}
fResiduals = new Double_t [fNPlanes];
+ fResidualsExclPlane = new Double_t [fNPlanes];
fWire_hit_did = new Double_t [fNPlanes];
fWire_hit_should = new Double_t [fNPlanes];
@@ -413,6 +414,7 @@ Int_t THcDC::DefineVariables( EMode mode )
{ "sp2_id", " (golden track) ", "fSp2_ID_best"},
{ "gtrack_nsp", " Number of space points in golden track ", "fNsp_best"},
{ "residual", "Residuals", "fResiduals"},
+ { "residualExclPlane", "Residuals", "fResidualsExclPlane"},
{ "wireHitDid","Wire did have matched track hit", "fWire_hit_did"},
{ "wireHitShould", "Wire should have matched track hit", "fWire_hit_should"},
{ 0 }
@@ -501,6 +503,7 @@ void THcDC::ClearEvent()
for(Int_t i=0;i<fNPlanes;i++) {
fResiduals[i] = 1000.0;
+ fResidualsExclPlane[i] = 1000.0;
fWire_hit_did[i] = 1000.0;
fWire_hit_should[i] = 1000.0;
}
@@ -646,6 +649,7 @@ void THcDC::SetFocalPlaneBestTrack(Int_t golden_track_index)
THcDCHit *hit = tr1->GetHit(ihit);
Int_t plane = hit->GetPlaneNum() - 1;
fResiduals[plane] = tr1->GetResidual(plane);
+ fResidualsExclPlane[plane] = tr1->GetResidualExclPlane(plane);
}
EfficiencyPerWire(golden_track_index);
}
@@ -1027,7 +1031,55 @@ void THcDC::TrackFit()
theDCTrack->SetVector(dray[0], dray[1], 0.0, dray[2], dray[3]);
}
theDCTrack->SetChisq(chi2);
+
+ // calculate ray without a plane in track
+ for(Int_t ipl_hit=0;ipl_hit < theDCTrack->GetNHits();ipl_hit++) {
+ if(theDCTrack->GetNFree() > 0) {
+ 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++) {
+ if (ihit != ipl_hit) {
+ 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++) {
+ if (ihit != ipl_hit) {
+ AA[irayp][jrayp] += fPlaneCoeffs[planes[ihit]][raycoeffmap[irayp]]*
+ fPlaneCoeffs[planes[ihit]][raycoeffmap[jrayp]]/
+ pow(fSigma[planes[ihit]],2);
+ }
+ }
+ }
+ }
+ }
+ //
+ // Solve 4x4 equations
+ // Should check that it is invertable
+ TVectorD dray(NUM_FPRAY);
+ AA.Invert();
+ dray = AA*TT;
+ Double_t coord=0.0;
+ for(Int_t ir=0;ir<NUM_FPRAY;ir++) {
+ coord += fPlaneCoeffs[planes[ipl_hit]][raycoeffmap[ir]]*dray[ir];
+ // cout << "ir = " << ir << ", dray[ir] = " << dray[ir] << endl;
+ }
+ Double_t residual = coords[ipl_hit] - coord;
+ theDCTrack->SetResidualExclPlane(planes[ipl_hit], residual);
+ }
+ }
}
+ //Calculate residual without plane
// Calculate residuals for each chamber if in single stub mode
// and there was a track found in each chamber
diff --git a/src/THcDC.h b/src/THcDC.h
index 9420f73..dd37637 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -122,6 +122,7 @@ protected:
Int_t fNSp; // Number of space points
Int_t fNsp_best; // Number of space points for gloden track
Double_t* fResiduals; //[fNPlanes] Array of residuals
+ Double_t* fResidualsExclPlane; //[fNPlanes] Array of residuals with plane excluded
Double_t* fWire_hit_did; //[fNPlanes]
Double_t* fWire_hit_should; //[fNPlanes]
diff --git a/src/THcDCTrack.cxx b/src/THcDCTrack.cxx
index 777bed2..e2f75c9 100644
--- a/src/THcDCTrack.cxx
+++ b/src/THcDCTrack.cxx
@@ -13,6 +13,7 @@ THcDCTrack::THcDCTrack(Int_t nplanes) : fnSP(0), fNHits(0)
fHits.clear();
fCoords.resize(nplanes);
fResiduals.resize(nplanes);
+ fResidualsExclPlane.resize(nplanes);
fDoubleResiduals.resize(nplanes);
}
diff --git a/src/THcDCTrack.h b/src/THcDCTrack.h
index c0cb14d..e899e74 100644
--- a/src/THcDCTrack.h
+++ b/src/THcDCTrack.h
@@ -37,6 +37,7 @@ public:
Double_t GetCoord(Int_t ip) const {return fCoords[ip];}
Double_t GetResidual(Int_t ip) const {return fResiduals[ip];}
Double_t GetResidual1() const {return fResiduals[0];}
+ Double_t GetResidualExclPlane(Int_t ip) const {return fResidualsExclPlane[ip];}
void GetRay(Double_t *ray) const {ray[0]=fX_fp; ray[1]=fY_fp; ray[2]=fXp_fp; ray[3]=fYp_fp;}
Double_t GetX() const {return fX_fp;}
Double_t GetY() const {return fY_fp;}
@@ -48,6 +49,7 @@ public:
void SetNFree(Int_t nfree) {fNfree = nfree;}
void SetCoord(Int_t ip, Double_t coord) {fCoords[ip] = coord;}
void SetResidual(Int_t ip, Double_t coord) {fResiduals[ip] = coord;}
+ void SetResidualExclPlane(Int_t ip, Double_t coord) {fResidualsExclPlane[ip] = coord;}
void SetDoubleResidual(Int_t ip, Double_t coord) {fDoubleResiduals[ip] = coord;}
void SetVector(Double_t x, Double_t y, Double_t z,
Double_t xp, Double_t yp) {fX_fp = x; fY_fp=y; fZ_fp=z; fXp_fp=xp; fYp_fp=yp;}
@@ -78,6 +80,7 @@ protected:
//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> fResidualsExclPlane; /* Residual on each plane */
std::vector<Double_t> fDoubleResiduals; /* Residual on each plane for single stub mode */
Double_t fX_fp, fY_fp, fZ_fp;
Double_t fXp_fp, fYp_fp;
--
GitLab