From e76a21c30fcab182715ddf75911bbacb76301b51 Mon Sep 17 00:00:00 2001
From: "Stephen A. Wood" <saw@jlab.org>
Date: Mon, 3 Jun 2013 17:11:14 -0400
Subject: [PATCH] Add LinkStubs method to THcDC. Comment out some debugging.
Compiles and doesn't crash.
---
src/THcDC.cxx | 178 +++++++++++++++++++++++++++++++++++
src/THcDC.h | 35 +++----
src/THcDriftChamber.cxx | 53 ++++++-----
src/THcDriftChamber.h | 4 +-
src/THcDriftChamberPlane.cxx | 2 +-
src/THcSpacePoint.h | 5 +
6 files changed, 232 insertions(+), 45 deletions(-)
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index 4b18439..884f3ed 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -284,10 +284,18 @@ Int_t THcDC::ReadDatabase( const TDatime& date )
{"dc_central_wire", fCentralWire, kDouble, fNPlanes},
{"dc_plane_time_zero", fPlaneTimeZero, kDouble, fNPlanes},
{"dc_sigma", fSigma, kDouble, fNPlanes},
+ {"single_stub",&fSingleStub, kInt},
+ {"ntracks_max_fp", &fNTracksMaxFP, kInt},
+ {"xt_track_criterion", &fXtTrCriterion, kDouble},
+ {"yt_track_criterion", &fYtTrCriterion, kDouble},
+ {"xpt_track_criterion", &fXptTrCriterion, kDouble},
+ {"ypt_track_criterion", &fYptTrCriterion, kDouble},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+ if(fNTracksMaxFP <= 0) fNTracksMaxFP = 10;
+ // if(fNTracksMaxFP > HNRACKS_MAX) fNTracksMaxFP = NHTRACKS_MAX;
cout << "Plane counts:";
for(Int_t i=0;i<fNPlanes;i++) {
cout << " " << fNWires[i];
@@ -444,6 +452,8 @@ Int_t THcDC::CoarseTrack( TClonesArray& /* tracks */ )
fChambers[i]->CorrectHitTimes();
fChambers[i]->LeftRight();
}
+ // Now link the stubs between chambers
+ LinkStubs();
ApplyCorrections();
@@ -464,6 +474,174 @@ Int_t THcDC::FineTrack( TClonesArray& tracks )
return 0;
}
+void THcDC::LinkStubs()
+{
+ // The logic is
+ // 0) Put all space points in a single list
+ // 1) loop over all space points as seeds isp1
+ // 2) Check if this space point is all ready in a track
+ // 3) loop over all succeeding space pointss isp2
+ // 4) check if there is a track-criterion match
+ // either add to existing track
+ // or if there is another point in same chamber
+ // make a copy containing isp2 rather than
+ // other point in same chamber
+ // 5) If hsingle_stub is set, make a track of all single
+ // stubs.
+
+ std::vector<THcSpacePoint*> fSp;
+ Int_t fNSp=0;
+ fSp.clear();
+ fSp.reserve(10);
+ // Make a vector of pointers to the SpacePoints
+ cout << "Linking " << fChambers[0]->GetNSpacePoints()
+ << " and " << fChambers[1]->GetNSpacePoints() << " stubs" << endl;
+ for(Int_t ich=0;ich<fNChambers;ich++) {
+ Int_t nchamber=fChambers[ich]->GetChamberNum();
+ TClonesArray* spacepointarray = fChambers[ich]->GetSpacePointsP();
+ for(Int_t isp=0;isp<fChambers[ich]->GetNSpacePoints();isp++) {
+ fSp.push_back(static_cast<THcSpacePoint*>(spacepointarray->At(isp)));
+ fSp[fNSp]->fNChamber = nchamber;
+ fNSp++;
+ }
+ }
+ Int_t ntracks_fp=0; // Number of Focal Plane tracks found
+ Double_t stubminx = 999999;
+ Double_t stubminy = 999999;
+ Double_t stubminxp = 999999;
+ Double_t stubminyp = 999999;
+ Int_t stub_tracks[MAXTRACKS];
+
+ if(!fSingleStub) {
+ for(Int_t isp1=0;isp1<fNSp-1;isp1++) {
+ Int_t sptracks=0;
+ // Now make sure this sp is not already used in a sp.
+ // Could this be done by having a sp point to the track it is in?
+ Int_t tryflag=1;
+ for(Int_t itrack=0;itrack<ntracks_fp;itrack++) {
+ for(Int_t isp=0;isp<fTrackSP[itrack].nSP;isp++) {
+ if(fTrackSP[itrack].spID[isp] == isp1) {
+ tryflag=0;
+ }
+ }
+ }
+ if(tryflag) { // SP not already part of a track
+ Int_t newtrack=1;
+ for(Int_t isp2=isp1+1;isp2<fNSp;isp2++) {
+ if(fSp[isp1]->fNChamber!=fSp[isp2]->fNChamber) {
+ Double_t *spstub1=fSp[isp1]->GetStubP();
+ Double_t *spstub2=fSp[isp2]->GetStubP();
+ Double_t dposx = spstub1[0] - spstub2[0];
+ Double_t dposy = spstub1[1] - spstub2[1];
+ Double_t dposxp = spstub1[2] - spstub2[2];
+ Double_t dposyp = spstub1[3] - spstub2[3];
+
+ // What is the point of saving these stubmin values. They
+ // Don't seem to be used anywhere except that they can be
+ // printed out if hbypass_track_eff_files is zero.
+ if(TMath::Abs(dposx)<TMath::Abs(stubminx)) stubminx = dposx;
+ if(TMath::Abs(dposy)<TMath::Abs(stubminy)) stubminy = dposy;
+ if(TMath::Abs(dposxp)<TMath::Abs(stubminxp)) stubminxp = dposxp;
+ if(TMath::Abs(dposyp)<TMath::Abs(stubminyp)) stubminyp = dposyp;
+
+ // if hbypass_track_eff_files == 0 then
+ // Print out each stubminX that is less that its criterion
+
+ if((TMath::Abs(dposx) < fXtTrCriterion)
+ && (TMath::Abs(dposy) < fYtTrCriterion)
+ && (TMath::Abs(dposxp) < fXptTrCriterion)
+ && (TMath::Abs(dposyp) < fYptTrCriterion)) {
+ if(newtrack) {
+ assert(sptracks==0);
+ //stubtest=1; Used in h_track_tests.f
+ // Make a new track if there are not to many
+ if(ntracks_fp < MAXTRACKS) {
+ sptracks=0; // Number of tracks with this seed
+ stub_tracks[sptracks++] = ntracks_fp;
+ fTrackSP[ntracks_fp].nSP=2;
+ fTrackSP[ntracks_fp].spID[0] = isp1;
+ fTrackSP[ntracks_fp].spID[1] = isp2;
+ // Now save the X, Y and XP for the two stubs
+ // in arrays hx_sp1, hy_sp1, hy_sp1, ... hxp_sp2
+ // Why not also YP?
+ // Skip for here. May be a diagnostic thing
+ newtrack = 0; // Make no more tracks in this loop
+ // (But could replace a SP?)
+ ntracks_fp++;
+ } else {
+ cout << "EPIC FAIL 1: Too many tracks found in THcDC::LinkStubs" << endl;
+ ntracks_fp=0;
+ // Do something here to fail this event
+ return;
+ }
+ } else {
+ // Check if there is another space point in the same chamber
+ for(Int_t itrack=0;itrack<sptracks;itrack++) {
+ Int_t track=stub_tracks[itrack];
+ Int_t spoint=0;
+ Int_t duppoint=0;
+ for(Int_t isp=0;fTrackSP[track].nSP;isp++) {
+ if(fSp[isp2]->fNChamber ==
+ fSp[fTrackSP[track].spID[isp]]->fNChamber) {
+ spoint=isp;
+ }
+ if(isp2==fTrackSP[track].spID[isp]) {
+ duppoint=1;
+ }
+ } // End loop over sp in tracks with isp1
+ // If there is no other space point in this chamber
+ // add this space point to current track(2)
+ if(!duppoint) {
+ if(!spoint) {
+ fTrackSP[track].spID[fTrackSP[track].nSP++] = isp2;
+ } else {
+ // If there is another point in the same chamber
+ // in this track create a new track with all the
+ // same space points except spoint
+ if(ntracks_fp < MAXTRACKS) {
+ stub_tracks[sptracks++] = ntracks_fp;
+ fTrackSP[ntracks_fp].nSP=fTrackSP[track].nSP;
+ for(Int_t isp=0;isp<fTrackSP[track].nSP;isp++) {
+ if(isp!=spoint) {
+ fTrackSP[ntracks_fp].spID[isp] = fTrackSP[track].spID[isp];
+ } else {
+ fTrackSP[ntracks_fp].spID[isp] = isp2;
+ } // End check for dup on copy
+ } // End copy of track
+ } else {
+ cout << "EPIC FAIL 2: Too many tracks found in THcDC::LinkStubs" << endl;
+ ntracks_fp=0;
+ // Do something here to fail this event
+ return; // Max # of allowed tracks
+ }
+ } // end if on same chamber
+ } // end if on duplicate point
+ } // end for over tracks with isp1
+ }
+ }
+ } // end test on same chamber
+ } // end isp2 loop over new space points
+ } // end test on tryflag
+ } // end isp1 outer loop over space points
+ } else { // Make track out of each single space point
+ for(Int_t isp=0;isp<fNSp;isp++) {
+ if(ntracks_fp<MAXTRACKS) {
+ fTrackSP[ntracks_fp].nSP=1;
+ fTrackSP[ntracks_fp].spID[0]=isp;
+ ntracks_fp++;
+ } else {
+ cout << "EPIC FAIL 3: Too many tracks found in THcDC::LinkStubs" << endl;
+ ntracks_fp=0;
+ // Do something here to fail this event
+ return; // Max # of allowed tracks
+ }
+ }
+ }
+ // Now list all hits on a track. What needs this
+ ///
+ ///
+ cout << "Found " << ntracks_fp << " tracks"<<endl;
+}
ClassImp(THcDC)
////////////////////////////////////////////////////////////////////////////////
diff --git a/src/THcDC.h b/src/THcDC.h
index e0b7b4d..12e9903 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -10,6 +10,7 @@
#include "THaTrackingDetector.h"
#include "THcHitList.h"
#include "THcRawDCHit.h"
+#include "THcSpacePoint.h"
#include "THcDriftChamberPlane.h"
#include "THcDriftChamber.h"
#include "TMath.h"
@@ -84,6 +85,12 @@ protected:
Double_t fNSperChan; /* TDC bin size */
Double_t fWireVelocity;
+ Int_t fSingleStub; /* If 1, single stubs make tracks */
+ Int_t fNTracksMaxFP;
+ Double_t fXtTrCriterion;
+ Double_t fYtTrCriterion;
+ Double_t fXptTrCriterion;
+ Double_t fYptTrCriterion;
// Each of these will be dimensioned with the number of chambers
Double_t* fXCenter;
@@ -113,32 +120,26 @@ protected:
Double_t* fPlaneTimeZero;
Double_t* fSigma;
+ // Useful derived quantities
+ // double tan_angle, sin_angle, cos_angle;
+
+ // Intermediate structure for building
+ static const char MAXTRACKS = 50;
+ struct TrackSP {
+ Int_t nSP; /* Number of space points in this track */
+ Int_t spID[10]; /* List of space points in this track */
+ } fTrackSP[MAXTRACKS];
+
std::vector<THcDriftChamberPlane*> fPlanes; // List of plane objects
std::vector<THcDriftChamber*> fChambers; // List of chamber objects
TClonesArray* fTrackProj; // projection of track onto scintillator plane
// and estimated match to TOF paddle
- // Useful derived quantities
- // double tan_angle, sin_angle, cos_angle;
-
- // static const char NDEST = 2;
- // struct DataDest {
- // Int_t* nthit;
- // Int_t* nahit;
- // Double_t* tdc;
- // Double_t* tdc_c;
- // Double_t* adc;
- // Double_t* adc_p;
- // Double_t* adc_c;
- // Double_t* offset;
- // Double_t* ped;
- // Double_t* gain;
- // } fDataDest[NDEST]; // Lookup table for decoder
-
void ClearEvent();
void DeleteArrays();
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
+ void LinkStubs();
void Setup(const char* name, const char* description);
diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx
index fe9b490..97c32bc 100644
--- a/src/THcDriftChamber.cxx
+++ b/src/THcDriftChamber.cxx
@@ -73,8 +73,8 @@ THcDriftChamber::THcDriftChamber( ) :
{
// Constructor
fPlanes.clear();
+ fSpacePoints = NULL;
}
-
//_____________________________________________________________________________
Int_t THcDriftChamber::Decode( const THaEvData& evdata )
{
@@ -271,50 +271,50 @@ Int_t THcDriftChamber::FindSpacePoints( void )
Int_t ip=thishit->GetPlaneNum(); // This is the absolute plane mumber
if(ip==YPlaneNum) yplane_hitind = ihit;
if(ip==YPlanePNum) yplanep_hitind = ihit;
- cout << " hit = " << ihit << " " << ip <<" " << thishit->GetWireNum()<<" " << thishit->GetPos()<< " " << thishit->GetPlaneNum()<< endl;
+ // cout << " hit = " << ihit << " " << ip <<" " << thishit->GetWireNum()<<" " << thishit->GetPos()<< " " << thishit->GetPlaneNum()<< endl;
}
// fPlanes is the array of planes for this chamber.
// cout << fPlanes[YPlaneInd]->GetNHits() << " "
// << fPlanes[YPlanePInd]->GetNHits() << " " << endl;
- cout << " Yplane ind " << YPlaneInd << " " << YPlanePInd << endl;
+ // cout << " Yplane ind " << YPlaneInd << " " << YPlanePInd << endl;
if (fPlanes[YPlaneInd]->GetNHits() == 1 && fPlanes[YPlanePInd]->GetNHits() == 1) {
- cout <<" Yplane info :" << fHits[yplane_hitind]->GetWireNum() << " "
- << fHits[yplane_hitind]->GetPos() << " "
- << fHits[yplanep_hitind]->GetWireNum() << " "
- << fHits[yplanep_hitind]->GetPos() << " "
- << fSpacePointCriterion << endl;
+ // cout <<" Yplane info :" << fHits[yplane_hitind]->GetWireNum() << " "
+ // << fHits[yplane_hitind]->GetPos() << " "
+ // << fHits[yplanep_hitind]->GetWireNum() << " "
+ // << fHits[yplanep_hitind]->GetPos() << " "
+ // << fSpacePointCriterion << endl;
}
if(fPlanes[YPlaneInd]->GetNHits() == 1 && fPlanes[YPlanePInd]->GetNHits() == 1
&& TMath::Abs(fHits[yplane_hitind]->GetPos() - fHits[yplanep_hitind]->GetPos())
< fSpacePointCriterion
&& fNhits <= 6) { // An easy case, probably one hit per plane
- cout << " call FindEasySpacePoint" << endl;
+ // cout << " call FindEasySpacePoint" << endl;
fEasySpacePoint = FindEasySpacePoint(yplane_hitind, yplanep_hitind);
}
if(!fEasySpacePoint) { // It's not easy
- cout << " call FindHardSpacePoints" << endl;
- cout << " nhits = " << fNhits << " " << fPlanes[YPlaneInd]->GetNHits() << " " << fPlanes[YPlanePInd]->GetNHits() << endl;
+ // cout << " call FindHardSpacePoints" << endl;
+ // cout << " nhits = " << fNhits << " " << fPlanes[YPlaneInd]->GetNHits() << " " << fPlanes[YPlanePInd]->GetNHits() << endl;
FindHardSpacePoints();
}
// We have our space points for this chamber
- cout << fNSpacePoints << " Space Points found" << endl;
+ // cout << fNSpacePoints << " Space Points found" << endl;
if(fNSpacePoints > 0) {
if(fRemove_Sppt_If_One_YPlane == 1) {
// The routine is specific to HMS
Int_t ndest=DestroyPoorSpacePoints();
- cout << ndest << " Poor space points destroyed" << endl;
+ // cout << ndest << " Poor space points destroyed" << endl;
// Loop over space points and remove those with less than 4 planes
// hit and missing hits in Y,Y' planes
}
- if(fNSpacePoints == 0) cout << "DestroyPoorSpacePoints() killed SP" << endl;
+ // if(fNSpacePoints == 0) cout << "DestroyPoorSpacePoints() killed SP" << endl;
Int_t nadded=SpacePointMultiWire();
if (nadded) cout << nadded << " Space Points added with SpacePointMultiWire()" << endl;
ChooseSingleHit();
SelectSpacePoints();
- if(fNSpacePoints == 0) cout << "SelectSpacePoints() killed SP" << endl;
+ // if(fNSpacePoints == 0) cout << "SelectSpacePoints() killed SP" << endl;
}
- cout << fNSpacePoints << " Space Points remain" << endl;
+ // cout << fNSpacePoints << " Space Points remain" << endl;
// Add these space points to the total list of space points for the
// the DC package. Do this in THcDC.cxx.
#if 0
@@ -343,7 +343,7 @@ Int_t THcDriftChamber::FindEasySpacePoint(Int_t yplane_hitind,Int_t yplanep_hiti
// a space point.
Int_t easy_space_point=0;
- cout << "Doing Find Easy Space Point" << endl;
+ // cout << "Doing Find Easy Space Point" << endl;
Double_t yt = (fHits[yplane_hitind]->GetPos() + fHits[yplanep_hitind]->GetPos())/2.0;
Double_t xt = 0.0;
Int_t num_xhits = 0;
@@ -353,7 +353,7 @@ Int_t THcDriftChamber::FindEasySpacePoint(Int_t yplane_hitind,Int_t yplanep_hiti
THcDCHit* thishit = fHits[ihit];
if(ihit!=yplane_hitind && ihit!=yplanep_hitind) { // x-like hit
// ysp and xsp are from h_generate_geometry
- cout << ihit << " " << thishit->GetPos() << " " << yt << " " << thishit->GetWirePlane()->GetYsp() << " " << thishit->GetWirePlane()->GetXsp() << endl;
+ // cout << ihit << " " << thishit->GetPos() << " " << yt << " " << thishit->GetWirePlane()->GetYsp() << " " << thishit->GetWirePlane()->GetXsp() << endl;
x_pos[ihit] = (thishit->GetPos()
-yt*thishit->GetWirePlane()->GetYsp())
/thishit->GetWirePlane()->GetXsp();
@@ -365,18 +365,18 @@ Int_t THcDriftChamber::FindEasySpacePoint(Int_t yplane_hitind,Int_t yplanep_hiti
}
Double_t max_dist = TMath::Sqrt(fSpacePointCriterion/2.);
xt = (num_xhits>0?xt/num_xhits:0.0);
- cout << " mean x = "<< xt << " max dist = " << max_dist << endl;
+ // cout << " mean x = "<< xt << " max dist = " << max_dist << endl;
easy_space_point = 1; // Assume we have an easy space point
// Rule it out if x points don't cluster well enough
for(Int_t ihit=0;ihit<fNhits;ihit++) {
- cout << "Hit " << ihit << " " << x_pos[ihit] << " " << TMath::Abs(xt-x_pos[ihit]) << endl;
+ // cout << "Hit " << ihit << " " << x_pos[ihit] << " " << TMath::Abs(xt-x_pos[ihit]) << endl;
if(ihit!=yplane_hitind && ihit!=yplanep_hitind) { // x-like hit
if(TMath::Abs(xt-x_pos[ihit]) >= max_dist)
{ easy_space_point=0; break;}
}
}
if(easy_space_point) { // Register the space point
- cout << "Easy Space Point " << xt << " " << yt << endl;
+ // cout << "Easy Space Point " << xt << " " << yt << endl;
THcSpacePoint* sp = (THcSpacePoint*)fSpacePoints->ConstructedAt(fNSpacePoints++);
sp->Clear();
sp->SetXY(xt, yt);
@@ -792,8 +792,8 @@ void THcDriftChamber::SelectSpacePoints()
}
}
}
- if(sp_count < fNSpacePoints)
- cout << "Reduced from " << fNSpacePoints << " to " << sp_count << " space points" << endl;
+ // if(sp_count < fNSpacePoints)
+ // cout << "Reduced from " << fNSpacePoints << " to " << sp_count << " space points" << endl;
fNSpacePoints = sp_count;
}
@@ -1049,13 +1049,13 @@ Double_t THcDriftChamber::FindStub(Int_t nhits, const std::vector<THcDCHit*>* hi
/fSigma[plane_list[ihit]];
}
}
- fAA3Inv[bitpat]->Print();
- cout << TT[0] << " " << TT[1] << " " << TT[2] << endl;
+ // fAA3Inv[bitpat]->Print();
+ // cout << TT[0] << " " << TT[1] << " " << TT[2] << endl;
// TT->Print();
//dstub = *(fAA3Inv[bitpat]) * TT;
TT *= *fAA3Inv[bitpat];
- cout << TT[0] << " " << TT[1] << " " << TT[2] << endl;
+ // cout << TT[0] << " " << TT[1] << " " << TT[2] << endl;
// cout << dstub[0] << " " << dstub[1] << " " << dstub[2] << endl;
// Calculate Chi2. Remember one power of sigma is in fStubCoefs
@@ -1079,6 +1079,7 @@ THcDriftChamber::~THcDriftChamber()
{
// Destructor. Remove variables from global list.
+ cout << fChamberNum << ": delete fSpacePoints: " << hex << fSpacePoints << endl;
delete fSpacePoints;
// Should delete the matrices
diff --git a/src/THcDriftChamber.h b/src/THcDriftChamber.h
index d50b3df..43c2dd3 100644
--- a/src/THcDriftChamber.h
+++ b/src/THcDriftChamber.h
@@ -47,10 +47,12 @@ public:
Int_t GetNSpacePoints() const { return(fNSpacePoints);}
Int_t GetNTracks() const { return fTrackProj->GetLast()+1; }
const TClonesArray* GetTrackHits() const { return fTrackProj; }
+ TClonesArray* GetSpacePointsP() const { return(fSpacePoints);}
+ Int_t GetChamberNum() const { return fChamberNum;}
// friend class THaScCalib;
- THcDriftChamber(); // for ROOT I/O
+ THcDriftChamber(); // for ROOT I/O // Why do we need this?
protected:
// Calibration
diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx
index 927e9cb..b68ab8a 100644
--- a/src/THcDriftChamberPlane.cxx
+++ b/src/THcDriftChamberPlane.cxx
@@ -107,6 +107,7 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
fParent = (THcDC*) GetParent();
// These are single variables here, but arrays in THcDriftChamber.
+ fSigma = fParent->GetSigma(fPlaneNum);
fChamberNum = fParent->GetNChamber(fPlaneNum);
fNWires = fParent->GetNWires(fPlaneNum);
fWireOrder = fParent->GetWireOrder(fPlaneNum);
@@ -118,7 +119,6 @@ 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();
diff --git a/src/THcSpacePoint.h b/src/THcSpacePoint.h
index 01e24fc..4122cb7 100644
--- a/src/THcSpacePoint.h
+++ b/src/THcSpacePoint.h
@@ -43,6 +43,10 @@ public:
void SetCombos(Int_t ncombos) { fNCombos=ncombos; };
Int_t GetCombos() { return fNCombos; };
+ // This is the chamber number (1,2), not index (0,1). Sometime
+ // we need figure out how to avoid confusion between number and index.
+ Int_t fNChamber;
+
protected:
Double_t fX;
@@ -51,6 +55,7 @@ public:
Int_t fNCombos;
std::vector<THcDCHit*> fHits;
Double_t fStub[4];
+ // Should we also have a pointer back to the chamber object
ClassDef(THcSpacePoint,0) // Drift Chamber class
};
--
GitLab