Skip to content
Snippets Groups Projects
Commit 294be162 authored by Stephen A. Wood's avatar Stephen A. Wood
Browse files

Implement more DC tracking. 

   Add track coordinates to track structure.
   Add code for single stubs
   Try to deal with plane number/index confusion
parent a8e4705a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/THcDC.cxx
+ 140
− 18
View file @ 294be162
...@@ -299,7 +299,7 @@ Int_t THcDC::ReadDatabase( const TDatime& date ) ...@@ -299,7 +299,7 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
{0} {0}
}; };
gHcParms->LoadParmValues((DBRequest*)&list,prefix); gHcParms->LoadParmValues((DBRequest*)&list,prefix);
fDebugDC=0; fDebugDC=1;
if(fNTracksMaxFP <= 0) fNTracksMaxFP = 10; if(fNTracksMaxFP <= 0) fNTracksMaxFP = 10;
// if(fNTracksMaxFP > HNRACKS_MAX) fNTracksMaxFP = NHTRACKS_MAX; // if(fNTracksMaxFP > HNRACKS_MAX) fNTracksMaxFP = NHTRACKS_MAX;
if (fDebugDC) cout << "Plane counts:"; if (fDebugDC) cout << "Plane counts:";
...@@ -461,6 +461,22 @@ Int_t THcDC::CoarseTrack( TClonesArray& /* tracks */ ) ...@@ -461,6 +461,22 @@ Int_t THcDC::CoarseTrack( TClonesArray& /* tracks */ )
// Now link the stubs between chambers // Now link the stubs between chambers
LinkStubs(); LinkStubs();
if(ntracks_fp > 0) TrackFit(); if(ntracks_fp > 0) TrackFit();
// Check for internal TrackFit errors
// Histogram the focal plane tracks
// Histograms made in h_fill_dc_fp_hist
// The following are one hist per track
// x_fp
// y_fp
// xp_fp
// yp_fp
// log chi2
// reduced chi2
// For each plane:
// double residual
// single residual
// Will need to make a track class that has all these things. Need to
// move the structure out of THcDC into it's own class which should probably
// inherit from a podd track class
ApplyCorrections(); ApplyCorrections();
...@@ -560,7 +576,7 @@ void THcDC::LinkStubs() ...@@ -560,7 +576,7 @@ void THcDC::LinkStubs()
&& (TMath::Abs(dposxp) < fXptTrCriterion) && (TMath::Abs(dposxp) < fXptTrCriterion)
&& (TMath::Abs(dposyp) < fYptTrCriterion)) { && (TMath::Abs(dposyp) < fYptTrCriterion)) {
if(newtrack) { if(newtrack) {
if (fDebugDC) cout << " new track" << endl; if (fDebugDC) cout << " new track" << endl;
assert(sptracks==0); assert(sptracks==0);
//stubtest=1; Used in h_track_tests.f //stubtest=1; Used in h_track_tests.f
// Make a new track if there are not to many // Make a new track if there are not to many
...@@ -584,13 +600,13 @@ void THcDC::LinkStubs() ...@@ -584,13 +600,13 @@ void THcDC::LinkStubs()
return; return;
} }
} else { } else {
if (fDebugDC) cout << " check if another space point in same chamber" << endl; if (fDebugDC) cout << " check if another space point in same chamber" << endl;
// Check if there is another space point in the same chamber // Check if there is another space point in the same chamber
for(Int_t itrack=0;itrack<sptracks;itrack++) { for(Int_t itrack=0;itrack<sptracks;itrack++) {
Int_t track=stub_tracks[itrack]; Int_t track=stub_tracks[itrack];
Int_t spoint=0; Int_t spoint=0;
Int_t duppoint=0; Int_t duppoint=0;
for(Int_t isp=0;fTrackSP[track].nSP;isp++) { for(Int_t isp=0;isp<fTrackSP[track].nSP;isp++) {
if(fSp[isp2]->fNChamber == if(fSp[isp2]->fNChamber ==
fSp[fTrackSP[track].spID[isp]]->fNChamber) { fSp[fTrackSP[track].spID[isp]]->fNChamber) {
spoint=isp; spoint=isp;
...@@ -690,6 +706,7 @@ void THcDC::TrackFit() ...@@ -690,6 +706,7 @@ void THcDC::TrackFit()
// Double_t chi2 = dummychi2; // Double_t chi2 = dummychi2;
// Int_t htrack_fit_num = itrack; // Int_t htrack_fit_num = itrack;
fTrackSP[itrack].fNfree = fTrackSP[itrack].fNHits - NUM_FPRAY; fTrackSP[itrack].fNfree = fTrackSP[itrack].fNHits - NUM_FPRAY;
Double_t chi2 = dummychi2;
if(fTrackSP[itrack].fNfree > 0) { if(fTrackSP[itrack].fNfree > 0) {
TVectorD TT(NUM_FPRAY); TVectorD TT(NUM_FPRAY);
TMatrixD AA(NUM_FPRAY,NUM_FPRAY); TMatrixD AA(NUM_FPRAY,NUM_FPRAY);
...@@ -697,7 +714,7 @@ void THcDC::TrackFit() ...@@ -697,7 +714,7 @@ void THcDC::TrackFit()
TT[irayp] = 0.0; TT[irayp] = 0.0;
for(Int_t ihit=0;ihit < fTrackSP[itrack].fNHits;ihit++) { for(Int_t ihit=0;ihit < fTrackSP[itrack].fNHits;ihit++) {
THcDCHit* hit=fTrackSP[itrack].fHits[ihit]; THcDCHit* hit=fTrackSP[itrack].fHits[ihit];
Int_t plane=hit->GetPlaneIndex(); Int_t plane=hit->GetPlaneNum()-1;
TT[irayp] += (hit->GetCoord()* TT[irayp] += (hit->GetCoord()*
fPlaneCoeffs[plane][raycoeffmap[irayp]]) fPlaneCoeffs[plane][raycoeffmap[irayp]])
/pow(fSigma[plane],2); /pow(fSigma[plane],2);
...@@ -711,7 +728,7 @@ void THcDC::TrackFit() ...@@ -711,7 +728,7 @@ void THcDC::TrackFit()
} else { } else {
for(Int_t ihit=0;ihit < fTrackSP[itrack].fNHits;ihit++) { for(Int_t ihit=0;ihit < fTrackSP[itrack].fNHits;ihit++) {
THcDCHit* hit=fTrackSP[itrack].fHits[ihit]; THcDCHit* hit=fTrackSP[itrack].fHits[ihit];
Int_t plane=hit->GetPlaneIndex(); Int_t plane=hit->GetPlaneNum()-1;
AA[irayp][jrayp] += fPlaneCoeffs[plane][raycoeffmap[irayp]]* AA[irayp][jrayp] += fPlaneCoeffs[plane][raycoeffmap[irayp]]*
fPlaneCoeffs[plane][raycoeffmap[jrayp]]/ fPlaneCoeffs[plane][raycoeffmap[jrayp]]/
pow(fSigma[plane],2); pow(fSigma[plane],2);
...@@ -729,33 +746,138 @@ void THcDC::TrackFit() ...@@ -729,33 +746,138 @@ void THcDC::TrackFit()
// if(bad_determinant) { // if(bad_determinant) {
// dray[0] = dray[1] = 10000.; dray[2] = dray[3] = 2.0; // dray[0] = dray[1] = 10000.; dray[2] = dray[3] = 2.0;
// } // }
Double_t chi2 = 0.0;
// Calculate hit coordinate for each plane for chi2 and efficiency // Calculate hit coordinate for each plane for chi2 and efficiency
// calculations // calculations
Double_t hdc_track_coord[fNPlanes]; fTrackSP[itrack].fCoords.clear();
fTrackSP[itrack].fResiduals.clear();
fTrackSP[itrack].fDoubleResiduals.clear();
for(Int_t iplane=0;iplane < fNPlanes; iplane++) { for(Int_t iplane=0;iplane < fNPlanes; iplane++) {
hdc_track_coord[iplane] = 0.0;// (itrk,pln) fTrackSP[itrack].fCoords.push_back(0.0);
fTrackSP[itrack].fResiduals.push_back(1000.0);
fTrackSP[itrack].fDoubleResiduals.push_back(1000.0);
for(Int_t ir=0;ir<NUM_FPRAY;ir++) { for(Int_t ir=0;ir<NUM_FPRAY;ir++) {
hdc_track_coord[iplane] += fPlaneCoeffs[iplane][raycoeffmap[ir]]*dray[ir]; fTrackSP[itrack].fCoords[iplane] += fPlaneCoeffs[iplane][raycoeffmap[ir]]*dray[ir];
} }
} }
if (fDebugDC) { // Compute Chi2 and residuals
cout << "Residuals:" << endl; chi2 = 0.0;
fTrackSP[itrack].fResiduals.clear();
for(Int_t ihit=0;ihit < fTrackSP[itrack].fNHits;ihit++) { for(Int_t ihit=0;ihit < fTrackSP[itrack].fNHits;ihit++) {
THcDCHit* hit=fTrackSP[itrack].fHits[ihit]; THcDCHit* hit=fTrackSP[itrack].fHits[ihit];
Int_t plane=hit->GetPlaneIndex(); Int_t plane=hit->GetPlaneNum()-1;
Double_t residual = hit->GetCoord() - hdc_track_coord[plane]; Double_t residual = hit->GetCoord() - fTrackSP[itrack].fCoords[plane];
cout << " " << plane << ": " << residual << endl; fTrackSP[itrack].fResiduals[plane] = residual;
chi2 += pow(residual/fSigma[plane],2); chi2 += pow(residual/fSigma[plane],2);
} }
if (fDebugDC) {
cout << "Residuals:" << endl;
for(Int_t ihit=0;ihit < fTrackSP[itrack].fNHits;ihit++) {
THcDCHit* hit=fTrackSP[itrack].fHits[ihit];
Int_t plane=hit->GetPlaneNum()-1;
cout << " " << plane << ": " << fTrackSP[itrack].fResiduals[ihit] << endl;
}
} }
fTrackSP[itrack].x_fp = dray[0];
fTrackSP[itrack].y_fp = dray[1];
fTrackSP[itrack].z_fp = 0.0;
fTrackSP[itrack].xp_fp = dray[2];
fTrackSP[itrack].xp_fp = dray[3];
} }
fTrackSP[itrack].chi2_fp = chi2;
} }
// Calculate residuals for each chamber if in single stub mode
// and there was a track found in each chamber
// Specific for two chambers. Can/should it be generalized?
if(fSingleStub != 0) {
if(ntracks_fp == 2) {
Int_t itrack=0;
Int_t ihit=0;
THcDCHit* hit=fTrackSP[itrack].fHits[ihit];
Int_t plane=hit->GetPlaneNum()-1;
Int_t chamber=fNChamber[plane];
if(chamber==1) {
itrack=1;
hit=fTrackSP[itrack].fHits[ihit];
plane=hit->GetPlaneNum()-1;
chamber=fNChamber[plane];
if(chamber==2) {
Double_t ray1[4];
Double_t ray2[4];
ray1[0] = fTrackSP[0].x_fp;
ray1[1] = fTrackSP[0].y_fp;
ray1[2] = fTrackSP[0].xp_fp;
ray1[3] = fTrackSP[0].yp_fp;
ray2[0] = fTrackSP[1].x_fp;
ray2[1] = fTrackSP[1].y_fp;
ray2[2] = fTrackSP[1].xp_fp;
ray2[3] = fTrackSP[1].yp_fp;
itrack = 1;
// Loop over hits in second chamber
for(Int_t ihit=0;ihit < fTrackSP[itrack].fNHits;ihit++) {
// Calculate residual in second chamber from first chamber track
THcDCHit* hit=fTrackSP[itrack].fHits[ihit];
Int_t plane=hit->GetPlaneNum()-1;
Double_t pos = DpsiFun(ray1,plane);
fTrackSP[itrack-1].fDoubleResiduals[plane] = hit->GetCoord() - pos;
// hdc_dbl_res(pln) = hdc_double_residual(1,pln) for hists
}
itrack=0;
// Loop over hits in first chamber
for(Int_t ihit=0;ihit < fTrackSP[itrack].fNHits;ihit++) {
// Calculate residual in first chamber from second chamber track
THcDCHit* hit=fTrackSP[itrack].fHits[ihit];
Int_t plane=hit->GetPlaneNum()-1;
Double_t pos = DpsiFun(ray1,plane);
fTrackSP[itrack+1].fDoubleResiduals[plane] = hit->GetCoord() - pos;
// hdc_dbl_res(pln) = hdc_double_residual(1,pln) for hists
}
}
}
}
}
// print tracks if hdebugtrackprint is on
} }
Double_t THcDC::DpsiFun(Double_t ray[4], Int_t plane)
{
/*
this function calculates the psi coordinate of the intersection
of a ray (defined by ray) with a hms wire chamber plane. the geometry
of the plane is contained in the coeff array calculated in the
array hplane_coeff
Note it is call by MINUIT via H_FCNCHISQ and so uses double precision
variables
the ray is defined by
x = (z-zt)*tan(xp) + xt
y = (z-zt)*tan(yp) + yt
at some fixed value of zt*
ray(1) = xt
ray(2) = yt
ray(3) = tan(xp)
ray(4) = tan(yp)
*/
Double_t infinity = 1.0E+20;
Double_t cinfinity = 1/infinity;
Double_t DpsiFun =
ray[2]*ray[1]*fPlaneCoeffs[plane][0] +
ray[3]*ray[0]*fPlaneCoeffs[plane][1] +
ray[2]*fPlaneCoeffs[plane][2] +
ray[3]*fPlaneCoeffs[plane][3] +
ray[0]*fPlaneCoeffs[plane][4] +
ray[1]*fPlaneCoeffs[plane][5];
Double_t denom = ray[2]*fPlaneCoeffs[plane][6]
+ ray[3]*fPlaneCoeffs[plane][7]
+ fPlaneCoeffs[plane][8];
if(TMath::Abs(denom) < cinfinity) {
DpsiFun = infinity;
} else {
DpsiFun = DpsiFun/denom;
}
return(DpsiFun);
}
ClassImp(THcDC) ClassImp(THcDC)
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
src/THcDC.h
+ 7
− 0
View file @ 294be162
...@@ -135,6 +135,12 @@ protected: ...@@ -135,6 +135,12 @@ protected:
Int_t fNHits; Int_t fNHits;
Int_t fNfree; /* Number of degrees of freedom */ Int_t fNfree; /* Number of degrees of freedom */
std::vector<THcDCHit*> 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 */
Double_t x_fp, y_fp, z_fp;
Double_t xp_fp, yp_fp;
Double_t chi2_fp;
} fTrackSP[MAXTRACKS]; } fTrackSP[MAXTRACKS];
std::vector<THcDriftChamberPlane*> fPlanes; // List of plane objects std::vector<THcDriftChamberPlane*> fPlanes; // List of plane objects
...@@ -148,6 +154,7 @@ protected: ...@@ -148,6 +154,7 @@ protected:
virtual Int_t DefineVariables( EMode mode = kDefine ); virtual Int_t DefineVariables( EMode mode = kDefine );
void LinkStubs(); void LinkStubs();
void TrackFit(); void TrackFit();
Double_t DpsiFun(Double_t ray[4], Int_t plane);
void Setup(const char* name, const char* description); void Setup(const char* name, const char* description);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment