diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index 26163a942b568ceab8363c9f9e9a61b07469e547..e79868cd26090cc462a36fdb37b40d3942064365 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -222,7 +222,7 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
fPitch = new Double_t [fNPlanes];
fCentralWire = new Double_t [fNPlanes];
fPlaneTimeZero = new Double_t [fNPlanes];
-
+ fSigma = new Double_t [fNPlanes];
DBRequest list[]={
{"dc_tdc_time_per_channel",&fNSperChan, kDouble},
@@ -250,6 +250,7 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
{"dc_pitch", fPitch, kDouble, fNPlanes},
{"dc_central_wire", fCentralWire, kDouble, fNPlanes},
{"dc_plane_time_zero", fPlaneTimeZero, kDouble, fNPlanes},
+ {"dc_sigma", fSigma, kDouble, fNPlanes},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
diff --git a/src/THcDC.h b/src/THcDC.h
index aac69980c7ea9708e8b71e14d59ae4aa339acd5f..7a69711388ef1c5685a7dd47174e4060e39868ce 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -44,6 +44,7 @@ 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 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];}
Double_t GetGammaAngle(Int_t plane) const { return fGammaAngle[plane-1];}
@@ -56,6 +57,7 @@ public:
Int_t GetDriftTimeSign(Int_t plane) const { return fDriftTimeSign[plane-1];}
Double_t GetPlaneTimeZero(Int_t plane) const { return fPlaneTimeZero[plane-1];}
+ Double_t GetSigma(Int_t plane) const { return fSigma[plane-1];}
Double_t GetNSperChan() const { return fNSperChan;}
@@ -109,6 +111,7 @@ protected:
Double_t* fPitch;
Double_t* fCentralWire;
Double_t* fPlaneTimeZero;
+ Double_t* fSigma;
THcDriftChamberPlane** fPlanes; // List of plane objects
THcDriftChamber** fChambers; // List of chamber objects
diff --git a/src/THcDCHit.h b/src/THcDCHit.h
index c44e53206e89c85f0332155bd441d6f08447e22f..e32cc6159bd9590170df39d4ebb15a175aebcec3 100644
--- a/src/THcDCHit.h
+++ b/src/THcDCHit.h
@@ -37,6 +37,7 @@ public:
Double_t GetTime() const { return fTime; }
Double_t GetDist() const { return fDist; }
Double_t GetPos() const { return fWire->GetPos(); } //Position of hit wire
+ Double_t GetCoord() const { return fCoord; }
Double_t GetdDist() const { return fdDist; }
Int_t GetCorrectedStatus() const { return fCorrected; }
@@ -47,6 +48,7 @@ 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 SetdDist(Double_t ddist) { fdDist = ddist; }
void SetFitDist(Double_t dist) { ftrDist = dist; }
Int_t GetPlaneNum() const { return fWirePlane->GetPlaneNum(); }
@@ -62,6 +64,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
+ Double_t fCoord; // Actual coordinate of hit
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?
diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx
index 7e66401b66981959460d8170c3b1c0141079f508..59312b75dd42c826c949ed2d2169f33224831529 100644
--- a/src/THcDriftChamber.cxx
+++ b/src/THcDriftChamber.cxx
@@ -19,6 +19,7 @@
#include "THaTrack.h"
#include "TClonesArray.h"
#include "TMath.h"
+#include "TVectorD.h"
#include "THaTrackProj.h"
@@ -26,6 +27,7 @@
#include <cstdio>
#include <cstdlib>
#include <iostream>
+#include <iomanip>
using namespace std;
@@ -79,7 +81,7 @@ Int_t THcDriftChamber::Decode( const THaEvData& evdata )
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcDriftChamber::Init( const TDatime& date )
{
- static const char* const here = "Init()";
+ // static const char* const here = "Init()";
Setup(GetName(), GetTitle());
@@ -120,6 +122,8 @@ Int_t THcDriftChamber::ReadDatabase( const TDatime& date )
DBRequest list[]={
{"_remove_sppt_if_one_y_plane",&fRemove_Sppt_If_One_YPlane, kInt},
{"dc_wire_velocity", &fWireVelocity, kDouble},
+ {"SmallAngleApprox", &fSmallAngleApprox, kInt},
+ {"stub_max_xpdiff", &fStubMaxXPDiff, kDouble},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
@@ -146,6 +150,64 @@ Int_t THcDriftChamber::ReadDatabase( const TDatime& date )
YPlanePNum = 11;
}
+ // Generate the HAA3INV matrix for all the acceptable combinations
+ // of hit planes. Try to make it as generic as possible
+ // pindex=0 -> Plane 1 missing, pindex5 -> plane 6 missing. Won't
+ // replicate the exact values used in the ENGINE, because the engine
+ // had one big list of matrices for both chambers, while here we will
+ // have a list just for one chamber. Also, call pindex, pmindex as
+ // we tend to use pindex as a plane index.
+ fCosBeta = new Double_t [fNPlanes];
+ fSinBeta = new Double_t [fNPlanes];
+ fTanBeta = new Double_t [fNPlanes];
+ fSigma = new Double_t [fNPlanes];
+ fPsi0 = new Double_t [fNPlanes];
+ fStubCoefs = new Double_t* [fNPlanes];
+ Int_t allplanes=0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ fCosBeta[ip] = TMath::Cos(fPlanes[ip]->GetBeta());
+ fSinBeta[ip] = TMath::Sin(fPlanes[ip]->GetBeta());
+ fTanBeta[ip] = fSinBeta[ip]/fCosBeta[ip];
+ fSigma[ip] = fPlanes[ip]->GetSigma();
+ fPsi0[ip] = fPlanes[ip]->GetPsi0();
+ fStubCoefs[ip] = fPlanes[ip]->GetStubCoef();
+ allplanes |= 1<<ip;
+ }
+ // Unordered map introduced in C++-11
+ // Can use unordered_map if using C++-11
+ // May not want to use map a all for performance, but using it now
+ // for code clarity
+ for(Int_t ipm1=0;ipm1<fNPlanes+1;ipm1++) { // Loop over missing plane1
+ for(Int_t ipm2=ipm1;ipm2<fNPlanes+1;ipm2++) {
+ if(ipm1==ipm2 && ipm1<fNPlanes) continue;
+ TMatrixDSym* AA3 = new TMatrixDSym(3);
+ for(Int_t i=0;i<3;i++) {
+ for(Int_t j=i;j<3;j++) {
+ (*AA3)[j][i] = 0.0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if(ipm1 != ip && ipm2 != ip) {
+ (*AA3)[j][i] += fStubCoefs[ip][i]*fStubCoefs[ip][j];
+ }
+ }
+ }
+ }
+ // Use the bit pattern of missing planes as an hash index.
+ // Perhaps it is more efficient to just use a regular array instead of map
+ Int_t bitpat = allplanes & ~(1<<ipm1) & ~(1<<ipm2);
+ // Should check that it is invertable
+ AA3->Invert();
+ fAA3Inv[bitpat] = AA3;
+ }
+ }
+
+#if 0
+ for(map<int,TMatrixDSym*>::iterator pm=fHaa3map.begin();
+ pm != fHaa3map.end(); pm++) {
+ cout << setbase(8) << pm->first << endl;
+ fAA3Inv[pm->first]->Print();
+ }
+#endif
+
fIsInit = true;
return kOK;
}
@@ -751,6 +813,235 @@ void THcDriftChamber::CorrectHitTimes()
}
}
}
+UInt_t THcDriftChamber::Count1Bits(UInt_t x)
+// From "Hacker's Delight"
+{
+ x = x - ((x >> 1) & 0x55555555);
+ x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
+ x = x + (x >> 8);
+ x = x + (x >> 16);
+ return x & 0x0000003F;
+}
+
+//_____________________________________________________________________________
+// HMS Specific
+void THcDriftChamber::LeftRight()
+{
+ // For each space point,
+ // Fit stubs to all possible left-right combinations of drift distances
+ // and choose the set with the minimum chi**2.
+
+ for(Int_t isp=0; isp<fNSpacePoints; isp++) {
+ // Build a bit pattern of which planes are hit
+ Int_t nhits = fSpacePoints[isp].nhits;
+ UInt_t bitpat = 0; // Bit pattern of which planes are hit
+ Double_t minchi2 = 1.0e10;
+ Double_t tmp_minchi2;
+ Double_t minxp = 0.25;
+ Int_t hasy1 = -1;
+ Int_t hasy2 = -1;
+ Int_t plusminusknown[nhits];
+ Int_t plusminusbest[nhits];
+ Int_t plusminus[nhits]; // ENGINE makes this array float. Why?
+ Int_t tmp_plusminus[nhits];
+ Int_t plane_list[nhits];
+ Double_t stub[4];
+ Double_t tmp_stub[4];
+
+ if(nhits < 0) {
+ cout << "THcDriftChamber::LeftRight() nhits < 0" << endl;
+ } else if (nhits==0) {
+ cout << "THcDriftChamber::LeftRight() nhits = 0" << endl;
+ }
+ for(Int_t ihit=0;ihit < nhits;ihit++) {
+ THcDCHit* hit = fSpacePoints[isp].hits[ihit];
+ Int_t pindex = hit->GetPlaneIndex();
+ plane_list[ihit] = pindex;
+
+ bitpat |= 1<<pindex;
+
+ plusminusknown[ihit] = 0;
+
+ if(pindex == YPlaneInd) hasy1 = ihit;
+ if(pindex == YPlanePInd) hasy2 = ihit;
+ }
+ Int_t smallAngOK = (hasy1>=0) && (hasy2>=0);
+ if(fSmallAngleApprox !=0 && smallAngOK) { // to small Angle L/R for Y,Y' planes
+ if(fSpacePoints[isp].hits[hasy1]->GetPos() <=
+ fSpacePoints[isp].hits[hasy2]->GetPos()) {
+ plusminusknown[hasy1] = -1;
+ plusminusknown[hasy2] = 1;
+ } else {
+ plusminusknown[hasy1] = 1;
+ plusminusknown[hasy2] = -1;
+ }
+ }
+ if(nhits < 2) {
+ cout << "THcDriftChamber::LeftRight: numhits-2 < 0" << endl;
+ } else if (nhits == 2) {
+ cout << "THcDriftChamber::LeftRight: numhits-2 = 0" << endl;
+ }
+ Int_t nplaneshit = Count1Bits(bitpat);
+ Int_t nplusminus = pow(2,nhits-2);
+
+ // Use bit value of integer word to set + or -
+ // Loop over all combinations of left right.
+ for(Int_t pmloop=0;pmloop<nplusminus;pmloop++) {
+ Int_t iswhit = 1;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ if(plusminusknown[ihit]!=0) {
+ plusminus[ihit] = plusminusknown[ihit];
+ } else {
+ // Max hits per point has to be less than 32.
+ if(pmloop & iswhit) {
+ plusminus[ihit] = 1;
+ } else {
+ plusminus[ihit] = -1;
+ }
+ iswhit <<= 1;
+ }
+ }
+ if (nplaneshit >= fNPlanes-1) {
+ Double_t chi2 = FindStub(nhits, fSpacePoints[isp].hits,
+ plane_list, bitpat, plusminus, stub);
+
+ //if(debugging)
+ //cout << "pmloop=" << pmloop << " Chi2=" << chi2 << endl;
+ // Take best chi2 IF x' of the stub agrees with x' as expected from x.
+ // Sometimes an incorrect x' gives a good chi2 for the stub, even though it is
+ // not the correct left/right combination for the real track.
+ // Rotate x'(=stub(3)) to hut coordinates and compare to x' expected from x.
+ // THIS ASSUMES STANDARD HMS TUNE!!!!, for which x' is approx. x/875.
+ if(chi2 < minchi2) {
+ if(stub[2]*fTanBeta[plane_list[0]]==-1.0) {
+ cout << "THcDriftChamber::LeftRight() Error 3" << endl;
+ }
+ Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]]
+ /(1+stub[2]*fTanBeta[plane_list[0]]);
+ // Tune depdendent. Definitely HMS specific
+ Double_t xp_expect = fSpacePoints[isp].x/875.0;
+ if(TMath::Abs(xp_fit-xp_expect)<fStubMaxXPDiff) {
+ minchi2 = chi2;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ plusminusbest[ihit] = plusminus[ihit];
+ }
+ for(Int_t i=0;i<4;i++) {
+ fSpacePoints[isp].stub[i] = stub[i];
+ }
+ } else { // Record best stub failing angle cut
+ tmp_minchi2 = chi2;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ tmp_plusminus[ihit] = plusminus[ihit];
+ }
+ for(Int_t i=0;i<4;i++) {
+ tmp_stub[i] = stub[i];
+ }
+ }
+ }
+ } else if (nplaneshit >= fNPlanes-2) { // Two planes missing
+ Double_t chi2 = FindStub(nhits, fSpacePoints[isp].hits,
+ plane_list, bitpat, plusminus, stub);
+ //if(debugging)
+ //cout << "pmloop=" << pmloop << " Chi2=" << chi2 << endl;
+ // Isn't this a bad idea, doing == with reals
+ if(stub[2]*fTanBeta[plane_list[0]] == -1.0) {
+ cout << "THcDriftChamber::LeftRight() Error 3" << endl;
+ }
+ Double_t xp_fit=stub[2]-fTanBeta[plane_list[0]]
+ /(1+stub[2]*fTanBeta[plane_list[0]]);
+ if(TMath::Abs(xp_fit) <= minxp) {
+ minxp = TMath::Abs(xp_fit);
+ minchi2 = chi2;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ plusminusbest[ihit] = plusminus[ihit];
+ }
+ for(Int_t i=0;i<4;i++) {
+ fSpacePoints[isp].stub[i] = stub[i];
+ }
+ }
+ } else {
+ cout << "Insufficient planes hit in THcDriftChamber::LeftRight()" << bitpat <<endl;
+ }
+ } // End loop of pm combinations
+
+ if(minchi2 > 9.9e9) { // No track passed angle cut
+ minchi2 = tmp_minchi2;
+ for(Int_t ihit=0;ihit<nhits;ihit++) {
+ plusminusbest[ihit] = tmp_plusminus[ihit];
+ }
+ for(Int_t i=0;i<4;i++) {
+ fSpacePoints[isp].stub[i] = tmp_stub[i];
+ }
+
+ }
+
+ // Calculate final coordinate based on plusminusbest
+ // Update the hit positions in the space points
+ for(Int_t ihit; ihit<nhits; ihit++) {
+ fSpacePoints[isp].hits[ihit]->SetLeftRight(plusminusbest[ihit]);
+ }
+
+ // Stubs are calculated in rotated coordinate system
+ // (I think this rotates in case chambers not perpendicular to central ray)
+ Int_t pindex=plane_list[0];
+ if(fSpacePoints[isp].stub[2] - fTanBeta[pindex] == -1.0) {
+ cout << "THcDriftChamber::LeftRight(): stub3 error" << endl;
+ }
+ stub[2] = (fSpacePoints[isp].stub[2] - fTanBeta[pindex])
+ / (1.0 + fSpacePoints[isp].stub[2]*fTanBeta[pindex]);
+ if(fSpacePoints[isp].stub[2]*fSinBeta[pindex] == -fCosBeta[pindex]) {
+ cout << "THcDriftChamber::LeftRight(): stub4 error" << endl;
+ }
+ stub[3] = fSpacePoints[isp].stub[3]
+ / (fSpacePoints[isp].stub[2]*fSinBeta[pindex]+fCosBeta[pindex]);
+ stub[0] = fSpacePoints[isp].stub[0]*fCosBeta[pindex]
+ - fSpacePoints[isp].stub[0]*stub[2]*fSinBeta[pindex];
+ stub[1] = fSpacePoints[isp].stub[1]
+ - fSpacePoints[isp].stub[1]*stub[3]*fSinBeta[pindex];
+ for(Int_t i=0;i<4;i++) {
+ fSpacePoints[isp].stub[i] = stub[i];
+ }
+ }
+ // Option to print stubs
+}
+ // if(fAA3Inv.find(bitpat) != fAAInv.end()) { // Valid hit combination
+//_____________________________________________________________________________
+Double_t THcDriftChamber::FindStub(Int_t nhits, THcDCHit** hits,
+ Int_t* plane_list, UInt_t bitpat,
+ Int_t* plusminus, Double_t* stub)
+{
+ // For a given combination of L/R, fit a stub to the space point
+ Double_t zeros[] = {0.0,0.0,0.0};
+ TVectorD TT; TT.Use(3, zeros);
+ TVectorD dstub;
+ Double_t dpos[3];
+
+ for(Int_t ihit=0;ihit<nhits; ihit++) {
+ dpos[ihit] = hits[ihit]->GetPos() + plusminus[ihit]*hits[ihit]->GetdDist()
+ - fPsi0[plane_list[ihit]];
+ for(Int_t index=0;index<3;index++) {
+ TT[index]+= dpos[ihit]*fStubCoefs[plane_list[ihit]][index]
+ /fSigma[plane_list[ihit]];
+ }
+ }
+ dstub = (*fAA3Inv[bitpat]) * TT;
+
+ // Calculate Chi2. Remember one power of sigma is in fStubCoefs
+ stub[0] = dstub[0];
+ stub[1] = dstub[1];
+ stub[2] = dstub[2];
+ stub[3] = 0.0;
+ Double_t chi2=0.0;
+ for(Int_t ihit=0;ihit<nhits; ihit++) {
+ chi2 += pow( dpos[ihit]/fSigma[plane_list[ihit]]
+ - fStubCoefs[plane_list[ihit]][0]*stub[0]
+ - fStubCoefs[plane_list[ihit]][1]*stub[1]
+ - fStubCoefs[plane_list[ihit]][2]*stub[2]
+ , 2);
+ }
+ return(chi2);
+}
+
//_____________________________________________________________________________
THcDriftChamber::~THcDriftChamber()
{
diff --git a/src/THcDriftChamber.h b/src/THcDriftChamber.h
index 0e7b6442d9a073aa44db9850c207f1a070de0651..db307dd06bfab3b204bf482e8c3820351d03da0b 100644
--- a/src/THcDriftChamber.h
+++ b/src/THcDriftChamber.h
@@ -10,6 +10,9 @@
#include "THaSubDetector.h"
#include "THcDriftChamberPlane.h"
#include "TClonesArray.h"
+#include "TMatrixDSym.h"
+
+#include <map>
#define MAX_SPACE_POINTS 50
#define MAX_HITS_PER_POINT 20
@@ -68,11 +71,19 @@ protected:
Int_t fMinCombos; // Minimum # pairs in a space point
Int_t fRemove_Sppt_If_One_YPlane;
Double_t fWireVelocity;
+ Int_t fSmallAngleApprox;
+ Double_t fStubMaxXPDiff;
Double_t fXCenter;
Double_t fYCenter;
Double_t fSpacePointCriterion;
Double_t fSpacePointCriterion2;
+ Double_t* fSinBeta;
+ Double_t* fCosBeta;
+ Double_t* fTanBeta;
+ Double_t* fSigma;
+ Double_t* fPsi0;
+ Double_t** fStubCoefs;
THcDriftChamberPlane* fPlanes[20]; // List of plane objects
@@ -90,6 +101,10 @@ protected:
Int_t SpacePointMultiWire(void);
void ChooseSingleHit(void);
void SelectSpacePoints(void);
+ UInt_t Count1Bits(UInt_t x);
+ void LeftRight(void);
+ Double_t FindStub(Int_t nhits, THcDCHit** hits, Int_t* plane_list,
+ UInt_t bitpat, Int_t* plusminus, Double_t* stub);
THcDCHit* fHits[10000]; /* All hits for this chamber */
// A simple structure until we figure out what we are doing.
@@ -99,11 +114,15 @@ protected:
Int_t nhits;
Int_t ncombos;
THcDCHit* hits[MAX_HITS_PER_POINT];
+ Double_t stub[4];
};
SpacePoint fSpacePoints[MAX_SPACE_POINTS];
Int_t fNSpacePoints;
Int_t fEasySpacePoint; /* This event is an easy space point */
+ Double_t* stubcoef[4];
+ std::map<int,TMatrixDSym*> fAA3Inv;
+
ClassDef(THcDriftChamber,0) // Drift Chamber class
};
diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx
index 8c63c5d9b3a03d385d6bc14c9938776c23e385d0..bd708bf6c51a7557ef3d09880390f5ab2c637586 100644
--- a/src/THcDriftChamberPlane.cxx
+++ b/src/THcDriftChamberPlane.cxx
@@ -112,14 +112,17 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
fCenter = fParent->GetCenter(fPlaneNum);
fCentralTime = fParent->GetCentralTime(fPlaneNum);
fDriftTimeSign = fParent->GetDriftTimeSign(fPlaneNum);
+ fSigma = fParent->GetSigma(fPlaneNum);
fNSperChan = fParent->GetNSperChan();
// Calculate Geometry Constants
// Do we want to move all this to the Chamber of DC Package leve
// as that is where these things will be needed?
+ Double_t z0 = fParent->GetZPos(fPlaneNum);
Double_t alpha = fParent->GetAlphaAngle(fPlaneNum);
Double_t beta = fParent->GetBetaAngle(fPlaneNum);
+ fBeta = beta;
Double_t gamma = fParent->GetGammaAngle(fPlaneNum);
Double_t cosalpha = TMath::Cos(alpha);
Double_t sinalpha = TMath::Sin(alpha);
@@ -127,13 +130,17 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
Double_t sinbeta = TMath::Sin(beta);
Double_t cosgamma = TMath::Cos(gamma);
Double_t singamma = TMath::Sin(gamma);
-
+
Double_t hzchi = -cosalpha*sinbeta + sinalpha*cosbeta*singamma;
Double_t hzpsi = sinalpha*sinbeta + cosalpha*cosbeta*singamma;
Double_t hxchi = -cosalpha*cosbeta - sinalpha*sinbeta*singamma;
Double_t hxpsi = sinalpha*cosbeta - cosalpha*sinbeta*singamma;
Double_t hychi = sinalpha*cosgamma;
Double_t hypsi = cosalpha*cosgamma;
+ Double_t stubxchi = -cosalpha;
+ Double_t stubxpsi = sinalpha;
+ Double_t stubychi = sinalpha;
+ Double_t stubypsi = cosalpha;
if(cosalpha <= 0.707) { // x-like wire, need dist from x=0 line
fReadoutX = 1;
@@ -146,9 +153,27 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
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;
- Double_t denom = sumsqupsi*sumsquchi-sumcross*sumcross;
- fXsp = hychi/denom;
- fYsp = -hxchi/denom;
+ Double_t denom1 = sumsqupsi*sumsquchi-sumcross*sumcross;
+ fPsi0 = (-z0*hzpsi*sumsquchi
+ +z0*hzchi*sumcross) / denom1;
+ Double_t hchi0 = (-z0*hzchi*sumsqupsi
+ +z0*hzpsi*sumcross) / denom1;
+ Double_t hphi0 = TMath::Sqrt(pow(z0+hzpsi*fPsi0+hzchi*hchi0,2)
+ + pow(hxpsi*fPsi0+hxchi*hchi0,2)
+ + pow(hypsi*fPsi0+hychi*hchi0,2) );
+ if(z0 < 0.0) hphi0 = -hphi0;
+
+ Double_t denom2 = stubxpsi*stubychi - stubxchi*stubypsi;
+
+ // Why are there 4, but only 3 used?
+ fStubCoef[0] = stubychi/(fSigma*denom2); // sin(a)/sigma
+ fStubCoef[1] = -stubxchi/(fSigma*denom2); // cos(a)/sigma
+ fStubCoef[2] = hphi0*fStubCoef[0]; // z0*sin(a)/sig
+ fStubCoef[3] = hphi0*fStubCoef[1]; // z0*cos(a)/sig
+
+ fXsp = hychi/denom2; // sin(a)
+ fYsp = -hxchi/denom2; // cos(a)
+
cout << fPlaneNum << " " << fNWires << " " << fWireOrder << endl;
diff --git a/src/THcDriftChamberPlane.h b/src/THcDriftChamberPlane.h
index e81c9ee5730e37a24169974bbea53fd405c0e0fc..475e9bc9d1d4038a4f1f8bede46e535f41d940b0 100644
--- a/src/THcDriftChamberPlane.h
+++ b/src/THcDriftChamberPlane.h
@@ -61,6 +61,10 @@ class THcDriftChamberPlane : public THaSubDetector {
Double_t GetReadoutCorr() { return fReadoutCorr; }
Double_t GetCentralTime() { return fCentralTime; }
Int_t GetDriftTimeSign() { return fDriftTimeSign; }
+ Double_t GetBeta() { return fBeta; }
+ Double_t GetSigma() { return fSigma; }
+ Double_t GetPsi0() { return fPsi0; }
+ Double_t* GetStubCoef() { return fStubCoef; }
protected:
@@ -81,6 +85,10 @@ class THcDriftChamberPlane : public THaSubDetector {
Double_t fPlaneTimeZero;
Double_t fXsp;
Double_t fYsp;
+ Double_t fSigma;
+ Double_t fPsi0;
+ Double_t fStubCoef[4];
+ Double_t fBeta;
Int_t fReadoutX;
Double_t fReadoutCorr;