diff --git a/Makefile b/Makefile
index fadfb2944800f96524642667e412b4977efa1e24..422e44b8456bc4d1847309026d5fb6da4dcf2d4b 100644
--- a/Makefile
+++ b/Makefile
@@ -15,7 +15,8 @@ SRC = src/THcInterface.cxx src/THcParmList.cxx src/THcAnalyzer.cxx \
src/THcSignalHit.cxx \
src/THcHodoscope.cxx src/THcScintillatorPlane.cxx \
src/THcHodoscopeHit.cxx \
- src/THcDriftChamber.cxx src/THcDriftChamberPlane.cxx \
+ src/THcDC.cxx src/THcDriftChamberPlane.cxx \
+ src/THcDriftChamber.cxx \
src/THcRawDCHit.cxx src/THcDCHit.cxx \
src/THcDCWire.cxx \
src/THcDCLookupTTDConv.cxx src/THcDCTimeToDistConv.cxx \
diff --git a/examples/hodtest.C b/examples/hodtest.C
index b92846d97e26933fda03fd850882f1a6bf5c9d6c..982e284adffdc261d98dec97cbf39ec5f105b017 100644
--- a/examples/hodtest.C
+++ b/examples/hodtest.C
@@ -38,7 +38,7 @@
// Add hodoscope
HMS->AddDetector( new THcHodoscope("hod", "Hodoscope" ));
HMS->AddDetector( new THcShower("cal", "Shower" ));
- HMS->AddDetector( new THcDriftChamber("dc", "Drift Chambers" ));
+ HMS->AddDetector( new THcDC("dc", "Drift Chambers" ));
THcAerogel* aerogel = new THcAerogel("aero", "Aerogel Cerenkov" );
HMS->AddDetector( aerogel );
diff --git a/examples/hodtest_cuts.def b/examples/hodtest_cuts.def
index 6ccbe52ceb55b33618ab5f59c576546ffe17aa64..f6003c2200a18b2f4d7156e51c0250c8b5ce2b06 100644
--- a/examples/hodtest_cuts.def
+++ b/examples/hodtest_cuts.def
@@ -9,6 +9,9 @@ RawDecode_master 1
Block: Decode
Decode_master !Pedestal_event
+Block: CoarseTracking
+CoarseTracking_master !Pedestal_event
+
Block: CoarseReconstruct
RawCoarseReconstruct !Pedestal_event
diff --git a/src/HallC_LinkDef.h b/src/HallC_LinkDef.h
index 608c0041beef3679cbdce6befbbf801ca61abd6d..ce520fa067dda33b10aebdad9bebadeb2336f719 100644
--- a/src/HallC_LinkDef.h
+++ b/src/HallC_LinkDef.h
@@ -18,6 +18,7 @@
#pragma link C++ class THcHodoscope+;
#pragma link C++ class THcScintillatorPlane+;
#pragma link C++ class THcHodoscopeHit+;
+#pragma link C++ class THcDC+;
#pragma link C++ class THcDriftChamber+;
#pragma link C++ class THcDriftChamberPlane+;
#pragma link C++ class THcRawDCHit+;
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..07eb1c6c3534741db4816b3df667b231e25018f9
--- /dev/null
+++ b/src/THcDC.cxx
@@ -0,0 +1,430 @@
+///////////////////////////////////////////////////////////////////////////////
+// //
+// THcDC //
+// //
+// Class for a generic hodoscope consisting of multiple //
+// planes with multiple paddles with phototubes on both ends. //
+// This differs from Hall A scintillator class in that it is the whole //
+// hodoscope array, not just one plane. //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THcDC.h"
+#include "THaEvData.h"
+#include "THaDetMap.h"
+#include "THcDetectorMap.h"
+#include "THcGlobals.h"
+#include "THcParmList.h"
+#include "VarDef.h"
+#include "VarType.h"
+#include "THaTrack.h"
+#include "TClonesArray.h"
+#include "TMath.h"
+
+#include "THaTrackProj.h"
+
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+
+using namespace std;
+
+//_____________________________________________________________________________
+THcDC::THcDC(
+ const char* name, const char* description,
+ THaApparatus* apparatus ) :
+ THaTrackingDetector(name,description,apparatus)
+{
+ // Constructor
+
+ // fTrackProj = new TClonesArray( "THaTrackProj", 5 );
+ fNPlanes = 0; // No planes until we make them
+
+}
+
+//_____________________________________________________________________________
+void THcDC::Setup(const char* name, const char* description)
+{
+
+ char prefix[2];
+ char parname[100];
+
+ THaApparatus *app = GetApparatus();
+ if(app) {
+ cout << app->GetName() << endl;
+ } else {
+ cout << "No apparatus found" << endl;
+ }
+
+ prefix[0]=tolower(app->GetName()[0]);
+ prefix[1]='\0';
+
+ string planenamelist;
+ DBRequest list[]={
+ {"dc_num_planes",&fNPlanes, kInt},
+ {"dc_num_chambers",&fNChambers, kInt},
+ {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
+ {"dc_wire_velocity",&fWireVelocity,kDouble},
+ {"dc_plane_names",&planenamelist, kString},
+ {0}
+ };
+
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+ cout << planenamelist << endl;
+ cout << "Drift Chambers: " << fNPlanes << " planes in " << fNChambers << " chambers" << endl;
+
+ vector<string> plane_names = vsplit(planenamelist);
+
+ if(plane_names.size() != (UInt_t) fNPlanes) {
+ cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names " << plane_names.size() << endl;
+ // Should quit. Is there an official way to quit?
+ }
+ fPlaneNames = new char* [fNPlanes];
+ for(Int_t i=0;i<fNPlanes;i++) {
+ fPlaneNames[i] = new char[plane_names[i].length()];
+ strcpy(fPlaneNames[i], plane_names[i].c_str());
+ }
+
+ char *desc = new char[strlen(description)+100];
+ fPlanes = new THcDriftChamberPlane* [fNPlanes];
+
+ for(Int_t i=0;i<fNPlanes;i++) {
+ strcpy(desc, description);
+ strcat(desc, " Plane ");
+ strcat(desc, fPlaneNames[i]);
+
+ fPlanes[i] = new THcDriftChamberPlane(fPlaneNames[i], desc, i+1, this);
+ cout << "Created Drift Chamber Plane " << fPlaneNames[i] << ", " << desc << endl;
+
+ }
+ fChambers = new THcDriftChamber* [fNChambers];
+ for(Int_t i=0;i<fNChambers;i++) {
+ sprintf(desc,"%s Chamber %d",description, i+1);
+
+ // Should construct a better chamber name
+ fChambers[i] = new THcDriftChamber(desc, desc, i+1, this);
+ cout << "Created Drift Chamber " << i+1 << ", " << desc << endl;
+
+
+ }
+}
+
+//_____________________________________________________________________________
+THcDC::THcDC( ) :
+ THaTrackingDetector()
+{
+ // Constructor
+}
+
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcDC::Init( const TDatime& date )
+{
+ static const char* const here = "Init()";
+
+ Setup(GetName(), GetTitle()); // Create the subdetectors here
+
+ // Should probably put this in ReadDatabase as we will know the
+ // maximum number of hits after setting up the detector map
+ THcHitList::InitHitList(fDetMap, "THcRawDCHit", 1000);
+
+ EStatus status;
+ // This triggers call of ReadDatabase and DefineVariables
+ if( (status = THaTrackingDetector::Init( date )) )
+ return fStatus=status;
+
+ // Initialize planes and add them to chambers
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if((status = fPlanes[ip]->Init( date ))) {
+ return fStatus=status;
+ } else {
+ Int_t chamber=fNChamber[ip];
+ fChambers[chamber-1]->AddPlane(fPlanes[ip]);
+ }
+ }
+ // Initialize chambers
+ for(Int_t ic=0;ic<fNChambers;ic++) {
+ if((status = fChambers[ic]->Init ( date ))) {
+ return fStatus=status;
+ }
+ }
+
+ // Replace with what we need for Hall C
+ // const DataDest tmp[NDEST] = {
+ // { &fRTNhit, &fRANhit, fRT, fRT_c, fRA, fRA_p, fRA_c, fROff, fRPed, fRGain },
+ // { &fLTNhit, &fLANhit, fLT, fLT_c, fLA, fLA_p, fLA_c, fLOff, fLPed, fLGain }
+ // };
+ // memcpy( fDataDest, tmp, NDEST*sizeof(DataDest) );
+
+ // Will need to determine which apparatus it belongs to and use the
+ // appropriate detector ID in the FillMap call
+ char EngineDID[4];
+
+ EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
+ EngineDID[1] = 'D';
+ EngineDID[2] = 'C';
+ EngineDID[3] = '\0';
+
+ if( gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0 ) {
+ Error( Here(here), "Error filling detectormap for %s.",
+ EngineDID);
+ return kInitError;
+ }
+
+ return fStatus = kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::ReadDatabase( const TDatime& date )
+{
+ // Read this detector's parameters from the database file 'fi'.
+ // This function is called by THaDetectorBase::Init() once at the
+ // beginning of the analysis.
+ // 'date' contains the date/time of the run being analyzed.
+
+ // static const char* const here = "ReadDatabase()";
+ char prefix[2];
+ char parname[100];
+
+ // Read data from database
+ // Pull values from the THcParmList instead of reading a database
+ // file like Hall A does.
+
+ // We will probably want to add some kind of method to gHcParms to allow
+ // bulk retrieval of parameters of interest.
+
+ // Will need to determine which spectrometer in order to construct
+ // the parameter names (e.g. hscin_1x_nr vs. sscin_1x_nr)
+
+ prefix[0]=tolower(GetApparatus()->GetName()[0]);
+
+ prefix[1]='\0';
+
+ fXCenter = new Double_t [fNChambers];
+ fYCenter = new Double_t [fNChambers];
+ fMinHits = new Int_t [fNChambers];
+ fMaxHits = new Int_t [fNChambers];
+ fMinCombos = new Int_t [fNChambers];
+ fSpace_Point_Criterion2 = new Double_t [fNChambers];
+
+ fTdcWinMin = new Int_t [fNPlanes];
+ fTdcWinMax = new Int_t [fNPlanes];
+ fCentralTime = new Int_t [fNPlanes];
+ fNWires = new Int_t [fNPlanes];
+ fNChamber = new Int_t [fNPlanes]; // Which chamber is this plane
+ fWireOrder = new Int_t [fNPlanes]; // Wire readout order
+ fDriftTimeSign = new Int_t [fNPlanes];
+
+ fZPos = new Double_t [fNPlanes];
+ fAlphaAngle = new Double_t [fNPlanes];
+ fBetaAngle = new Double_t [fNPlanes];
+ fGammaAngle = new Double_t [fNPlanes];
+ fPitch = new Double_t [fNPlanes];
+ fCentralWire = new Double_t [fNPlanes];
+ fPlaneTimeZero = new Double_t [fNPlanes];
+
+
+ DBRequest list[]={
+ {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
+ {"dc_wire_velocity",&fWireVelocity,kDouble},
+
+ {"dc_xcenter", fXCenter, kDouble, fNChambers},
+ {"dc_ycenter", fYCenter, kDouble, fNChambers},
+ {"min_hit", fMinHits, kInt, fNChambers},
+ {"max_pr_hits", fMaxHits, kInt, fNChambers},
+ {"min_combos", fMinCombos, kInt, fNChambers},
+ {"space_point_criterion", fSpace_Point_Criterion2, kDouble, fNChambers},
+
+ {"dc_tdc_min_win", fTdcWinMin, kInt, fNPlanes},
+ {"dc_tdc_max_win", fTdcWinMax, kInt, fNPlanes},
+ {"dc_central_time", fCentralTime, kInt, fNPlanes},
+ {"dc_nrwire", fNWires, kInt, fNPlanes},
+ {"dc_chamber_planes", fNChamber, kInt, fNPlanes},
+ {"dc_wire_counting", fWireOrder, kInt, fNPlanes},
+ {"dc_drifttime_sign", fDriftTimeSign, kInt, fNPlanes},
+
+ {"dc_zpos", fZPos, kDouble, fNPlanes},
+ {"dc_alpha_angle", fAlphaAngle, kDouble, fNPlanes},
+ {"dc_beta_angle", fBetaAngle, kDouble, fNPlanes},
+ {"dc_gamma_angle", fGammaAngle, kDouble, fNPlanes},
+ {"dc_pitch", fPitch, kDouble, fNPlanes},
+ {"dc_central_wire", fCentralWire, kDouble, fNPlanes},
+ {"dc_plane_time_zero", fPlaneTimeZero, kDouble, fNPlanes},
+ {0}
+ };
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+
+ cout << "Plane counts:";
+ for(Int_t i=0;i<fNPlanes;i++) {
+ cout << " " << fNWires[i];
+ }
+ cout << endl;
+
+ fIsInit = true;
+
+ return kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::DefineVariables( EMode mode )
+{
+ // Initialize global variables and lookup table for decoder
+
+ if( mode == kDefine && fIsSetup ) return kOK;
+ fIsSetup = ( mode == kDefine );
+
+ // Register variables in global list
+
+ RVarDef vars[] = {
+ { "nhit", "Number of DC hits", "fNhits" },
+ { 0 }
+ };
+ return DefineVarsFromList( vars, mode );
+
+}
+
+//_____________________________________________________________________________
+THcDC::~THcDC()
+{
+ // Destructor. Remove variables from global list.
+
+ if( fIsSetup )
+ RemoveVariables();
+ if( fIsInit )
+ DeleteArrays();
+ if (fTrackProj) {
+ fTrackProj->Clear();
+ delete fTrackProj; fTrackProj = 0;
+ }
+}
+
+//_____________________________________________________________________________
+void THcDC::DeleteArrays()
+{
+ // Delete member arrays. Used by destructor.
+
+ delete [] fNWires; fNWires = NULL;
+ // delete [] fSpacing; fSpacing = NULL;
+ // delete [] fCenter; fCenter = NULL; // This 2D. What is correct way to delete?
+
+ // delete [] fRA_c; fRA_c = NULL;
+ // delete [] fRA_p; fRA_p = NULL;
+ // delete [] fRA; fRA = NULL;
+ // delete [] fLA_c; fLA_c = NULL;
+ // delete [] fLA_p; fLA_p = NULL;
+ // delete [] fLA; fLA = NULL;
+ // delete [] fRT_c; fRT_c = NULL;
+ // delete [] fRT; fRT = NULL;
+ // delete [] fLT_c; fLT_c = NULL;
+ // delete [] fLT; fLT = NULL;
+
+ // delete [] fRGain; fRGain = NULL;
+ // delete [] fLGain; fLGain = NULL;
+ // delete [] fRPed; fRPed = NULL;
+ // delete [] fLPed; fLPed = NULL;
+ // delete [] fROff; fROff = NULL;
+ // delete [] fLOff; fLOff = NULL;
+ // delete [] fTWalkPar; fTWalkPar = NULL;
+ // delete [] fTrigOff; fTrigOff = NULL;
+
+ // delete [] fHitPad; fHitPad = NULL;
+ // delete [] fTime; fTime = NULL;
+ // delete [] fdTime; fdTime = NULL;
+ // delete [] fYt; fYt = NULL;
+ // delete [] fYa; fYa = NULL;
+}
+
+//_____________________________________________________________________________
+inline
+void THcDC::ClearEvent()
+{
+ // Reset per-event data.
+ fNhits = 0;
+
+ for(Int_t i=0;i<fNChambers;i++) {
+ fChambers[i]->Clear();
+ }
+
+
+ // fTrackProj->Clear();
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::Decode( const THaEvData& evdata )
+{
+
+ ClearEvent();
+
+ // Get the Hall C style hitlist (fRawHitList) for this event
+ fNhits = THcHitList::DecodeToHitList(evdata);
+
+ // Let each plane get its hits
+ Int_t nexthit = 0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
+ }
+
+ // Let each chamber get its hits
+ for(Int_t ic=0;ic<fNChambers;ic++) {
+ fChambers[ic]->ProcessHits();
+ }
+#if 0
+ // fRawHitList is TClones array of THcRawDCHit objects
+ for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
+ THcRawDCHit* hit = (THcRawDCHit *) fRawHitList->At(ihit);
+ // cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
+ // << endl;
+ for(Int_t imhit = 0; imhit < hit->fNHits; imhit++) {
+ // cout << " " << imhit << " " << hit->fTDC[imhit]
+ // << endl;
+ }
+ }
+ // cout << endl;
+#endif
+
+ return fNhits;
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::ApplyCorrections( void )
+{
+ return(0);
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::CoarseTrack( TClonesArray& /* tracks */ )
+{
+ // Calculation of coordinates of particle track cross point with scint
+ // plane in the detector coordinate system. For this, parameters of track
+ // reconstructed in THaVDC::CoarseTrack() are used.
+ //
+ // Apply corrections and reconstruct the complete hits.
+ //
+ // static const Double_t sqrt2 = TMath::Sqrt(2.);
+ for(Int_t i=0;i<fNChambers;i++) {
+
+ fChambers[i]->FindSpacePoints();
+ }
+
+ ApplyCorrections();
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcDC::FineTrack( TClonesArray& tracks )
+{
+ // Reconstruct coordinates of particle track cross point with scintillator
+ // plane, and copy the data into the following local data structure:
+ //
+ // Units of measurements are meters.
+
+ // Calculation of coordinates of particle track cross point with scint
+ // plane in the detector coordinate system. For this, parameters of track
+ // reconstructed in THaVDC::FineTrack() are used.
+
+ return 0;
+}
+
+ClassImp(THcDC)
+////////////////////////////////////////////////////////////////////////////////
diff --git a/src/THcDC.h b/src/THcDC.h
new file mode 100644
index 0000000000000000000000000000000000000000..d73491b49572d829720744d6f038eda2db06e7f0
--- /dev/null
+++ b/src/THcDC.h
@@ -0,0 +1,145 @@
+#ifndef ROOT_THcDC
+#define ROOT_THcDC
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// THcDC //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THaTrackingDetector.h"
+#include "THcHitList.h"
+#include "THcRawDCHit.h"
+#include "THcDriftChamberPlane.h"
+#include "THcDriftChamber.h"
+#include "TMath.h"
+
+//class THaScCalib;
+class TClonesArray;
+
+class THcDC : public THaTrackingDetector, public THcHitList {
+
+public:
+ THcDC( const char* name, const char* description = "",
+ THaApparatus* a = NULL );
+ virtual ~THcDC();
+
+ virtual Int_t Decode( const THaEvData& );
+ virtual EStatus Init( const TDatime& run_time );
+ virtual Int_t CoarseTrack( TClonesArray& tracks );
+ virtual Int_t FineTrack( TClonesArray& tracks );
+
+ virtual Int_t ApplyCorrections( void );
+
+ // Int_t GetNHits() const { return fNhit; }
+
+ Int_t GetNTracks() const { return fTrackProj->GetLast()+1; }
+ const TClonesArray* GetTrackHits() const { return fTrackProj; }
+
+ Int_t GetNWires(Int_t plane) const { return fNWires[plane-1];}
+ Int_t GetNChamber(Int_t plane) const { return fNChamber[plane-1];}
+ Int_t GetWireOrder(Int_t plane) const { return fWireOrder[plane-1];}
+ Int_t GetPitch(Int_t plane) const { return fPitch[plane-1];}
+ Int_t GetCentralWire(Int_t plane) const { return fCentralWire[plane-1];}
+ Int_t GetTdcWinMin(Int_t plane) const { return fTdcWinMin[plane-1];}
+ Int_t GetTdcWinMax(Int_t plane) const { return fTdcWinMax[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];}
+
+ Int_t GetMinHits(Int_t chamber) const { return fMinHits[chamber-1];}
+ Int_t GetMaxHits(Int_t chamber) const { return fMaxHits[chamber-1];}
+ Int_t GetMinCombos(Int_t chamber) const { return fMinCombos[chamber-1];}
+ Double_t GetSpacePointCriterion(Int_t chamber) const { return TMath::Sqrt(fSpace_Point_Criterion2[chamber-1]);}
+
+ Double_t GetPlaneTimeZero(Int_t plane) const { return fPlaneTimeZero[plane-1];}
+
+ Double_t GetNSperChan() const { return fNSperChan;}
+
+ Double_t GetCenter(Int_t plane) const {
+ Int_t chamber = GetNChamber(plane)-1;
+ return
+ fXCenter[chamber]*sin(fAlphaAngle[plane-1]) +
+ fYCenter[chamber]*cos(fAlphaAngle[plane-1]);
+ }
+ // friend class THaScCalib;
+
+ THcDC(); // for ROOT I/O
+protected:
+
+ // Calibration
+
+ // Per-event data
+ Int_t fNhits;
+
+ // Potential Hall C parameters. Mostly here for demonstration
+ Int_t fNPlanes;
+ char** fPlaneNames;
+ Int_t fNChambers;
+
+ Double_t fNSperChan; /* TDC bin size */
+ Double_t fWireVelocity;
+
+ // Each of these will be dimensioned with the number of chambers
+ Double_t* fXCenter;
+ Double_t* fYCenter;
+ Int_t* fMinHits;
+ Int_t* fMaxHits;
+ Int_t* fMinCombos;
+ Double_t* fSpace_Point_Criterion2;
+
+ // Each of these will be dimensioned with the number of planes
+ // A THcDCPlane class object will need to access the value for
+ // its plane number. Should we have a Get method for each or
+ Int_t* fTdcWinMin;
+ Int_t* fTdcWinMax;
+ Int_t* fCentralTime;
+ Int_t* fNWires; // Number of wires per plane
+ Int_t* fNChamber;
+ Int_t* fWireOrder;
+ Int_t* fDriftTimeSign;
+
+ Double_t* fZPos;
+ Double_t* fAlphaAngle;
+ Double_t* fBetaAngle;
+ Double_t* fGammaAngle;
+ Double_t* fPitch;
+ Double_t* fCentralWire;
+ Double_t* fPlaneTimeZero;
+
+ THcDriftChamberPlane** fPlanes; // List of plane objects
+ 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 Setup(const char* name, const char* description);
+
+ ClassDef(THcDC,0) // Drift Chamber class
+};
+
+////////////////////////////////////////////////////////////////////////////////
+
+#endif
diff --git a/src/THcDCHit.cxx b/src/THcDCHit.cxx
index 0484b1ad39d1f2d1ee0171dfe0d2f5194f6a0ce0..3b1523c444ce0d4088ab5b422359d7e3d1cfd9ac 100644
--- a/src/THcDCHit.cxx
+++ b/src/THcDCHit.cxx
@@ -9,10 +9,31 @@
#include "THcDCHit.h"
#include "THcDCTimeToDistConv.h"
+#include <iostream>
+
+using std::cout;
+using std::endl;
+
ClassImp(THcDCHit)
const Double_t THcDCHit::kBig = 1.e38; // Arbitrary large value
+//_____________________________________________________________________________
+void THcDCHit::Print( Option_t* opt ) const
+{
+ // Print hit info
+
+ cout << "Hit: wire=" << GetWireNum()
+ << "/" << (fWirePlane ? fWirePlane->GetName() : "??")
+ << " wpos=" << GetPos()
+ << " time=" << GetTime()
+ << " drift=" << GetDist();
+ // << " res=" << GetResolution()
+ // << " z=" << GetZ()
+ if( *opt != 'C' )
+ cout << endl;
+}
+
//_____________________________________________________________________________
Double_t THcDCHit::ConvertTimeToDist()
{
@@ -51,15 +72,17 @@ Int_t THcDCHit::Compare( const TObject* obj ) const
Int_t myWireNum = fWire->GetNum();
Int_t hitWireNum = hit->GetWire()->GetNum();
- // Compare wire numbers
+ Int_t myPlaneNum = GetPlaneNum();
+ Int_t hitPlaneNum = hit->GetPlaneNum();
+ if (myPlaneNum < hitPlaneNum) return -1;
+ if (myPlaneNum > hitPlaneNum) return 1;
+ // If planes are the same, compare wire numbers
if (myWireNum < hitWireNum) return -1;
if (myWireNum > hitWireNum) return 1;
- if (myWireNum == hitWireNum) {
- // If wire numbers are the same, compare times
- Double_t hitTime = hit->GetTime();
- if (fTime < hitTime) return -1;
- if (fTime > hitTime) return 1;
- }
+ // If wire numbers are the same, compare times
+ Double_t hitTime = hit->GetTime();
+ if (fTime < hitTime) return -1;
+ if (fTime > hitTime) return 1;
return 0;
}
diff --git a/src/THcDCHit.h b/src/THcDCHit.h
index 79f88e7123f7a8d39d805719e7ce2b74ec831316..fb886f6b1138ecea0d8ef53b757a4a15786c7d21 100644
--- a/src/THcDCHit.h
+++ b/src/THcDCHit.h
@@ -9,13 +9,17 @@
#include "TObject.h"
#include "THcDCWire.h"
+#include "THcDriftChamberPlane.h"
+#include "THcDriftChamber.h"
#include <cstdio>
class THcDCHit : public TObject {
public:
- THcDCHit( THcDCWire* wire=NULL, Int_t rawtime=0, Double_t time=0.0 ) :
- fWire(wire), fRawTime(rawtime), fTime(time), fDist(0.0), ftrDist(kBig) {
+ THcDCHit( THcDCWire* wire=NULL, Int_t rawtime=0, Double_t time=0.0,
+ THcDriftChamberPlane* wp=0) :
+ fWire(wire), fRawTime(rawtime), fTime(time), fWirePlane(wp),
+ fDist(0.0), ftrDist(kBig) {
ConvertTimeToDist();
}
virtual ~THcDCHit() {}
@@ -23,6 +27,7 @@ public:
virtual Double_t ConvertTimeToDist();
Int_t Compare ( const TObject* obj ) const;
Bool_t IsSortable () const { return kTRUE; }
+ virtual void Print( Option_t* opt="" ) const;
// Get and Set Functions
THcDCWire* GetWire() const { return fWire; }
@@ -33,22 +38,32 @@ public:
Double_t GetPos() const { return fWire->GetPos(); } //Position of hit wire
Double_t GetdDist() const { return fdDist; }
+ THcDriftChamberPlane* GetWirePlane() const { return fWirePlane; }
+
+
void SetWire(THcDCWire * wire) { fWire = wire; }
void SetRawTime(Int_t time) { fRawTime = time; }
void SetTime(Double_t time) { fTime = time; }
void SetDist(Double_t dist) { fDist = dist; }
void SetdDist(Double_t ddist) { fdDist = ddist; }
void SetFitDist(Double_t dist) { ftrDist = dist; }
+ Int_t GetPlaneNum() const { return fWirePlane->GetPlaneNum(); }
+ Int_t GetPlaneIndex() const { return fWirePlane->GetPlaneIndex(); }
+ Int_t GetChamberNum() const { return fWirePlane->GetChamberNum(); }
protected:
static const Double_t kBig; //!
- THcDCWire* fWire; // Wire on which the hit occurred
+ THcDCWire* fWire; // Wire on which the hit occurred
Int_t fRawTime; // TDC value (channels)
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 fdDist; // uncertainty in fDist (for chi2 calc)
Double_t ftrDist; // (Perpendicular) distance from the track
+
+ THcDriftChamber* fChamber; //! Pointer to parent wire plane
+
private:
THcDCHit( const THcDCHit& );
diff --git a/src/THcDriftChamber.cxx b/src/THcDriftChamber.cxx
index 1d53b20c58f52e7465596ea9240a58bd1755eb6b..1a77d4d7f230f59aaed0cc5349e44412130cb4dc 100644
--- a/src/THcDriftChamber.cxx
+++ b/src/THcDriftChamber.cxx
@@ -1,18 +1,17 @@
///////////////////////////////////////////////////////////////////////////////
// //
-// THcDriftChamber //
+// THcDriftChamber //
// //
-// Class for a generic hodoscope consisting of multiple //
-// planes with multiple paddles with phototubes on both ends. //
-// This differs from Hall A scintillator class in that it is the whole //
-// hodoscope array, not just one plane. //
+// Subdetector class to hold a bunch of planes constituting a chamber //
+// This class will be created by the THcDC class which will also create //
+// the plane objects. //
+// The THcDC class will then pass this class a list of the planes. //
// //
///////////////////////////////////////////////////////////////////////////////
#include "THcDriftChamber.h"
-#include "THaEvData.h"
-#include "THaDetMap.h"
-#include "THcDetectorMap.h"
+#include "THcDC.h"
+#include "THcDCHit.h"
#include "THcGlobals.h"
#include "THcParmList.h"
#include "VarDef.h"
@@ -33,14 +32,16 @@ using namespace std;
//_____________________________________________________________________________
THcDriftChamber::THcDriftChamber(
const char* name, const char* description,
- THaApparatus* apparatus ) :
- THaTrackingDetector(name,description,apparatus)
+ const Int_t chambernum, THaDetectorBase* parent ) :
+ THaSubDetector(name,description,parent)
{
// Constructor
// fTrackProj = new TClonesArray( "THaTrackProj", 5 );
fNPlanes = 0; // No planes until we make them
+ fChamberNum = chambernum;
+
}
//_____________________________________________________________________________
@@ -48,7 +49,6 @@ void THcDriftChamber::Setup(const char* name, const char* description)
{
char prefix[2];
- char parname[100];
THaApparatus *app = GetApparatus();
if(app) {
@@ -60,181 +60,92 @@ void THcDriftChamber::Setup(const char* name, const char* description)
prefix[0]=tolower(app->GetName()[0]);
prefix[1]='\0';
- string planenamelist;
- DBRequest list[]={
- {"dc_num_planes",&fNPlanes, kInt},
- {"dc_num_chambers",&fNChambers, kInt},
- {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
- {"dc_wire_velocity",&fWireVelocity,kDouble},
- {"dc_plane_names",&planenamelist, kString},
- {0}
- };
-
- gHcParms->LoadParmValues((DBRequest*)&list,prefix);
- cout << planenamelist << endl;
- cout << "Drift Chambers: " << fNPlanes << " planes in " << fNChambers << " chambers" << endl;
-
- vector<string> plane_names = vsplit(planenamelist);
-
- if(plane_names.size() != (UInt_t) fNPlanes) {
- cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names " << plane_names.size() << endl;
- // Should quit. Is there an official way to quit?
- }
- fPlaneNames = new char* [fNPlanes];
- for(Int_t i=0;i<fNPlanes;i++) {
- fPlaneNames[i] = new char[plane_names[i].length()];
- strcpy(fPlaneNames[i], plane_names[i].c_str());
- }
-
- char *desc = new char[strlen(description)+100];
- fPlanes = new THcDriftChamberPlane* [fNPlanes];
-
- for(Int_t i=0;i<fNPlanes;i++) {
- strcpy(desc, description);
- strcat(desc, " Plane ");
- strcat(desc, fPlaneNames[i]);
-
- fPlanes[i] = new THcDriftChamberPlane(fPlaneNames[i], desc, i+1, this);
- cout << "Created Drift Chamber Plane " << fPlaneNames[i] << ", " << desc << endl;
-
- }
-
}
//_____________________________________________________________________________
THcDriftChamber::THcDriftChamber( ) :
- THaTrackingDetector()
+ THaSubDetector()
{
// Constructor
}
+//_____________________________________________________________________________
+Int_t THcDriftChamber::Decode( const THaEvData& evdata )
+{
+ cout << "THcDriftChamber::Decode called" << endl;
+ return 0;
+}
+
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcDriftChamber::Init( const TDatime& date )
{
static const char* const here = "Init()";
- Setup(GetName(), GetTitle()); // Create the subdetectors here
+ Setup(GetName(), GetTitle());
- // Should probably put this in ReadDatabase as we will know the
- // maximum number of hits after setting up the detector map
- THcHitList::InitHitList(fDetMap, "THcRawDCHit", 1000);
-
EStatus status;
// This triggers call of ReadDatabase and DefineVariables
- if( (status = THaTrackingDetector::Init( date )) )
+ if( (status = THaSubDetector::Init( date )) )
return fStatus=status;
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- if((status = fPlanes[ip]->Init( date ))) {
- return fStatus=status;
- }
- }
-
- // Replace with what we need for Hall C
- // const DataDest tmp[NDEST] = {
- // { &fRTNhit, &fRANhit, fRT, fRT_c, fRA, fRA_p, fRA_c, fROff, fRPed, fRGain },
- // { &fLTNhit, &fLANhit, fLT, fLT_c, fLA, fLA_p, fLA_c, fLOff, fLPed, fLGain }
- // };
- // memcpy( fDataDest, tmp, NDEST*sizeof(DataDest) );
-
- // Will need to determine which apparatus it belongs to and use the
- // appropriate detector ID in the FillMap call
- char EngineDID[4];
+ return fStatus = kOK;
+}
- EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
- EngineDID[1] = 'D';
- EngineDID[2] = 'C';
- EngineDID[3] = '\0';
-
- if( gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0 ) {
- Error( Here(here), "Error filling detectormap for %s.",
- EngineDID);
- return kInitError;
+void THcDriftChamber::AddPlane(THcDriftChamberPlane *plane)
+{
+ cout << "Added plane " << plane->GetPlaneNum() << " to chamber " << fChamberNum << endl;
+ plane->SetPlaneIndex(fNPlanes);
+ fPlanes[fNPlanes] = plane;
+
+ // HMS Specific
+ // Check if this is a Y Plane
+ if(plane->GetPlaneNum() == YPlaneNum) {
+ YPlaneInd = fNPlanes;
+ } else if (plane->GetPlaneNum() == YPlanePNum) {
+ YPlanePInd = fNPlanes;
}
+ fNPlanes++;
- return fStatus = kOK;
+ return;
}
//_____________________________________________________________________________
Int_t THcDriftChamber::ReadDatabase( const TDatime& date )
{
- // Read this detector's parameters from the database file 'fi'.
- // This function is called by THaDetectorBase::Init() once at the
- // beginning of the analysis.
- // 'date' contains the date/time of the run being analyzed.
- // static const char* const here = "ReadDatabase()";
+ cout << "THcDriftChamber::ReadDatabase()" << endl;
char prefix[2];
- char parname[100];
-
- // Read data from database
- // Pull values from the THcParmList instead of reading a database
- // file like Hall A does.
-
- // We will probably want to add some kind of method to gHcParms to allow
- // bulk retrieval of parameters of interest.
-
- // Will need to determine which spectrometer in order to construct
- // the parameter names (e.g. hscin_1x_nr vs. sscin_1x_nr)
-
prefix[0]=tolower(GetApparatus()->GetName()[0]);
-
prefix[1]='\0';
-
- fXCenter = new Double_t [fNChambers];
- fYCenter = new Double_t [fNChambers];
-
- fTdcWinMin = new Int_t [fNPlanes];
- fTdcWinMax = new Int_t [fNPlanes];
- fCentralTime = new Int_t [fNPlanes];
- fNWires = new Int_t [fNPlanes];
- fNChamber = new Int_t [fNPlanes]; // Which chamber is this plane
- fWireOrder = new Int_t [fNPlanes]; // Wire readout order
- fDriftTimeSign = new Int_t [fNPlanes];
-
- fZPos = new Double_t [fNPlanes];
- fAlphaAngle = new Double_t [fNPlanes];
- fBetaAngle = new Double_t [fNPlanes];
- fGammaAngle = new Double_t [fNPlanes];
- fPitch = new Double_t [fNPlanes];
- fCentralWire = new Double_t [fNPlanes];
- fPlaneTimeZero = new Double_t [fNPlanes];
-
-
DBRequest list[]={
- {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
- {"dc_wire_velocity",&fWireVelocity,kDouble},
-
- {"dc_xcenter", fXCenter, kDouble, fNChambers},
- {"dc_ycenter", fYCenter, kDouble, fNChambers},
-
- {"dc_tdc_min_win", fTdcWinMin, kInt, fNPlanes},
- {"dc_tdc_max_win", fTdcWinMax, kInt, fNPlanes},
- {"dc_central_time", fCentralTime, kInt, fNPlanes},
- {"dc_nrwire", fNWires, kInt, fNPlanes},
- {"dc_chamber_planes", fNChamber, kInt, fNPlanes},
- {"dc_wire_counting", fWireOrder, kInt, fNPlanes},
- {"dc_drifttime_sign", fDriftTimeSign, kInt, fNPlanes},
-
- {"dc_zpos", fZPos, kDouble, fNPlanes},
- {"dc_alpha_angle", fAlphaAngle, kDouble, fNPlanes},
- {"dc_beta_angle", fBetaAngle, kDouble, fNPlanes},
- {"dc_gamma_angle", fGammaAngle, kDouble, fNPlanes},
- {"dc_pitch", fPitch, kDouble, fNPlanes},
- {"dc_central_wire", fCentralWire, kDouble, fNPlanes},
- {"dc_plane_time_zero", fPlaneTimeZero, kDouble, fNPlanes},
+ {"_remove_sppt_if_one_y_plane",&fRemove_Sppt_If_One_YPlane, kInt},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
- cout << "Plane counts:";
- for(Int_t i=0;i<fNPlanes;i++) {
- cout << " " << fNWires[i];
- }
- cout << endl;
+ // Get parameters parent knows about
+ THcDC* fParent;
+ fParent = (THcDC*) GetParent();
+ fMinHits = fParent->GetMinHits(fChamberNum);
+ fMaxHits = fParent->GetMaxHits(fChamberNum);
+ fMinCombos = fParent->GetMinCombos(fChamberNum);
- fIsInit = true;
+ cout << "Chamber " << fChamberNum << " Min/Max: " << fMinHits << "/" << fMaxHits << endl;
+ fSpacePointCriterion = fParent->GetSpacePointCriterion(fChamberNum);
+ fSpacePointCriterion2 = fSpacePointCriterion*fSpacePointCriterion;
+
+ // HMS Specific
+ // Hard code Y plane numbers. Should be able to get from wire angle
+ if(fChamberNum == 1) {
+ YPlaneNum = 2;
+ YPlanePNum = 5;
+ } else {
+ YPlaneNum = 8;
+ YPlanePNum = 11;
+ }
+
+ fIsInit = true;
return kOK;
}
@@ -249,165 +160,607 @@ Int_t THcDriftChamber::DefineVariables( EMode mode )
// Register variables in global list
// RVarDef vars[] = {
- // { "nlthit", "Number of Left paddles TDC times", "fLTNhit" },
- // { "nrthit", "Number of Right paddles TDC times", "fRTNhit" },
- // { "nlahit", "Number of Left paddles ADCs amps", "fLANhit" },
- // { "nrahit", "Number of Right paddles ADCs amps", "fRANhit" },
- // { "lt", "TDC values left side", "fLT" },
- // { "lt_c", "Corrected times left side", "fLT_c" },
- // { "rt", "TDC values right side", "fRT" },
- // { "rt_c", "Corrected times right side", "fRT_c" },
- // { "la", "ADC values left side", "fLA" },
- // { "la_p", "Corrected ADC values left side", "fLA_p" },
- // { "la_c", "Corrected ADC values left side", "fLA_c" },
- // { "ra", "ADC values right side", "fRA" },
- // { "ra_p", "Corrected ADC values right side", "fRA_p" },
- // { "ra_c", "Corrected ADC values right side", "fRA_c" },
- // { "nthit", "Number of paddles with l&r TDCs", "fNhit" },
- // { "t_pads", "Paddles with l&r coincidence TDCs", "fHitPad" },
- // { "y_t", "y-position from timing (m)", "fYt" },
- // { "y_adc", "y-position from amplitudes (m)", "fYa" },
- // { "time", "Time of hit at plane (s)", "fTime" },
- // { "dtime", "Est. uncertainty of time (s)", "fdTime" },
- // { "dedx", "dEdX-like deposited in paddle", "fAmpl" },
- // { "troff", "Trigger offset for paddles", "fTrigOff"},
- // { "trn", "Number of tracks for hits", "GetNTracks()" },
- // { "trx", "x-position of track in det plane", "fTrackProj.THaTrackProj.fX" },
- // { "try", "y-position of track in det plane", "fTrackProj.THaTrackProj.fY" },
- // { "trpath", "TRCS pathlen of track to det plane","fTrackProj.THaTrackProj.fPathl" },
- // { "trdx", "track deviation in x-position (m)", "fTrackProj.THaTrackProj.fdX" },
- // { "trpad", "paddle-hit associated with track", "fTrackProj.THaTrackProj.fChannel" },
+ // { "nhit", "Number of DC hits", "fNhits" },
// { 0 }
// };
// return DefineVarsFromList( vars, mode );
return kOK;
-}
-//_____________________________________________________________________________
-THcDriftChamber::~THcDriftChamber()
+}
+void THcDriftChamber::ProcessHits( void)
{
- // Destructor. Remove variables from global list.
+ // Make a list of hits for whole chamber
+ fNhits = 0;
- if( fIsSetup )
- RemoveVariables();
- if( fIsInit )
- DeleteArrays();
- if (fTrackProj) {
- fTrackProj->Clear();
- delete fTrackProj; fTrackProj = 0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ TClonesArray* hitsarray = fPlanes[ip]->GetHits();
+ for(Int_t ihit=0;ihit<fPlanes[ip]->GetNHits();ihit++) {
+ fHits[fNhits++] = static_cast<THcDCHit*>(hitsarray->At(ihit));
+ }
}
+ // cout << "ThcDriftChamber::ProcessHits() " << fNhits << " hits" << endl;
}
-//_____________________________________________________________________________
-void THcDriftChamber::DeleteArrays()
+
+Int_t THcDriftChamber::FindSpacePoints( void )
{
- // Delete member arrays. Used by destructor.
+ // Following h_pattern_recognition.f, but just for one chamber
- delete [] fNWires; fNWires = NULL;
- // delete [] fSpacing; fSpacing = NULL;
- // delete [] fCenter; fCenter = NULL; // This 2D. What is correct way to delete?
-
- // delete [] fRA_c; fRA_c = NULL;
- // delete [] fRA_p; fRA_p = NULL;
- // delete [] fRA; fRA = NULL;
- // delete [] fLA_c; fLA_c = NULL;
- // delete [] fLA_p; fLA_p = NULL;
- // delete [] fLA; fLA = NULL;
- // delete [] fRT_c; fRT_c = NULL;
- // delete [] fRT; fRT = NULL;
- // delete [] fLT_c; fLT_c = NULL;
- // delete [] fLT; fLT = NULL;
-
- // delete [] fRGain; fRGain = NULL;
- // delete [] fLGain; fLGain = NULL;
- // delete [] fRPed; fRPed = NULL;
- // delete [] fLPed; fLPed = NULL;
- // delete [] fROff; fROff = NULL;
- // delete [] fLOff; fLOff = NULL;
- // delete [] fTWalkPar; fTWalkPar = NULL;
- // delete [] fTrigOff; fTrigOff = NULL;
-
- // delete [] fHitPad; fHitPad = NULL;
- // delete [] fTime; fTime = NULL;
- // delete [] fdTime; fdTime = NULL;
- // delete [] fYt; fYt = NULL;
- // delete [] fYa; fYa = NULL;
+ // Code below is specifically for HMS chambers with Y and Y' planes
+
+ Int_t yplane_hitind=0;
+ Int_t yplanep_hitind=0;
+
+ fNSpacePoints=0;
+ fEasySpacePoint = 0;
+ // Should really build up array of hits
+ if(fNhits >= fMinHits && fNhits < fMaxHits) {
+ /* Has this list of hits already been cut the way it should have at this point? */
+
+ for(Int_t ihit=0;ihit<fNhits;ihit++) {
+ THcDCHit* thishit = fHits[ihit];
+ Int_t ip=thishit->GetPlaneNum(); // This is the absolute plane mumber
+ if(ip==YPlaneNum) yplane_hitind = ihit;
+ if(ip==YPlanePNum) yplanep_hitind = ihit;
+ }
+ // fPlanes is the array of planes for this chamber.
+ // cout << fPlanes[YPlaneInd]->GetNHits() << " "
+ // << fPlanes[YPlanePInd]->GetNHits() << " " << endl;
+ if (fPlanes[YPlaneInd]->GetNHits() == 1 && fPlanes[YPlanePInd]->GetNHits() == 1) {
+ cout << 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
+ fEasySpacePoint = FindEasySpacePoint(yplane_hitind, yplanep_hitind);
+ }
+ if(!fEasySpacePoint) { // It's not easy
+ FindHardSpacePoints();
+ }
+
+ // We have our space points for this chamber
+ 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;
+ // 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;
+ Int_t nadded=SpacePointMultiWire();
+ if (nadded) cout << nadded << " Space Points added with SpacePointMultiWire()" << endl;
+ ChooseSingleHit();
+ SelectSpacePoints();
+ if(fNSpacePoints == 0) cout << "SelectSpacePoints() killed SP" << 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
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ Int_t np = fPlanes[ip]->GetPlaneNum(); // Actuall plane number of this plane
+ // (0-11) or (1-12)?
+ TClonesArray* fHits = fPlanes[ip]->GetHits();
+
+ for(Int_t ihit=0;ihit<fNhits;ihit++) { // Looping over all hits in all planes of the chamber
+ THcDCHit* hit = static_cast<THcDCHit*>(fHits->At(ihit));
+ //
+
+ }
+ }
+#endif
+ }
+ return(fNSpacePoints);
}
//_____________________________________________________________________________
-inline
-void THcDriftChamber::ClearEvent()
+// HMS Specific
+Int_t THcDriftChamber::FindEasySpacePoint(Int_t yplane_hitind,Int_t yplanep_hitind)
{
- // Reset per-event data.
+ // Simplified HMS find_space_point routing. It is given all y hits and
+ // checks to see if all x-like hits are close enough together to make
+ // a space point.
+
+ Int_t easy_space_point=0;
+ 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;
+ Double_t x_pos[MAX_HITS_PER_POINT];
+
+ for(Int_t ihit=0;ihit<fNhits;ihit++) {
+ THcDCHit* thishit = fHits[ihit];
+ if(ihit!=yplane_hitind && ihit!=yplanep_hitind) { // x-like hit
+ // ysp and xsp are from h_generate_geometry
+ x_pos[ihit] = (thishit->GetPos()
+ -yt*thishit->GetWirePlane()->GetYsp())
+ /thishit->GetWirePlane()->GetXsp();
+ xt += x_pos[ihit];
+ num_xhits++;
+ } else {
+ x_pos[ihit] = 0.0;
+ }
+ }
+ xt = (num_xhits>0?xt/num_xhits:0.0);
+ cout << "xt=" << xt << 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] << " " << xt << endl;
+ if(ihit!=yplane_hitind && ihit!=yplanep_hitind) { // x-like hit
+ if(TMath::Abs(xt-x_pos[ihit]) >= fSpacePointCriterion)
+ { easy_space_point=0; break;}
+ }
+ }
+ if(easy_space_point) { // Register the space point
+ cout << "Easy Space Point " << xt << " " << yt << endl;
+ fNSpacePoints = 1;
+ fSpacePoints[0].x = xt;
+ fSpacePoints[0].y = yt;
+ fSpacePoints[0].nhits = fNhits;
+ fSpacePoints[0].ncombos = 0; // No combos
+ for(Int_t ihit=0;ihit<fNhits;ihit++) {
+ THcDCHit* thishit = fHits[ihit];
+ fSpacePoints[0].hits[ihit] = thishit;
+ }
+ }
+ return(easy_space_point);
+}
- fTrackProj->Clear();
+//_____________________________________________________________________________
+// Generic
+Int_t THcDriftChamber::FindHardSpacePoints()
+{
+ Int_t MAX_NUMBER_PAIRS=1000; // Where does this get set?
+ struct Pair {
+ THcDCHit* hit1;
+ THcDCHit* hit2;
+ Double_t x, y;
+ };
+ Pair pairs[MAX_NUMBER_PAIRS];
+ //
+ Int_t ntest_points=0;
+ for(Int_t ihit1=0;ihit1<fNhits-1;ihit1++) {
+ THcDCHit* hit1=fHits[ihit1];
+ THcDriftChamberPlane* plane1 = hit1->GetWirePlane();
+ for(Int_t ihit2=ihit1+1;ihit2<fNhits;ihit2++) {
+ if(ntest_points < MAX_NUMBER_PAIRS) {
+ THcDCHit* hit2=fHits[ihit2];
+ THcDriftChamberPlane* plane2 = hit2->GetWirePlane();
+ Double_t determinate = plane1->GetXsp()*plane2->GetYsp()
+ -plane1->GetYsp()*plane2->GetXsp();
+ if(TMath::Abs(determinate) > 0.3) { // 0.3 is sin(alpha1-alpha2)=sin(17.5)
+ pairs[ntest_points].hit1 = hit1;
+ pairs[ntest_points].hit2 = hit2;
+ pairs[ntest_points].x = (hit1->GetPos()*plane2->GetYsp()
+ - hit2->GetPos()*plane1->GetYsp())
+ /determinate;
+ pairs[ntest_points].y = (hit2->GetPos()*plane1->GetXsp()
+ - hit1->GetPos()*plane2->GetXsp())
+ /determinate;
+ ntest_points++;
+ }
+ }
+ }
+ }
+ Int_t ncombos=0;
+ struct Combo {
+ Pair* pair1;
+ Pair* pair2;
+ };
+ Combo combos[10*MAX_NUMBER_PAIRS];
+ for(Int_t ipair1=0;ipair1<ntest_points-1;ipair1++) {
+ for(Int_t ipair2=ipair1+1;ipair2<ntest_points;ipair2++) {
+ Double_t dist2 = pow(pairs[ipair1].x - pairs[ipair2].x,2)
+ + pow(pairs[ipair1].y - pairs[ipair2].y,2);
+ if(dist2 <= fSpacePointCriterion2) {
+ combos[ncombos].pair1 = &pairs[ipair1];
+ combos[ncombos].pair2 = &pairs[ipair2];
+ ncombos++;
+ }
+ }
+ }
+ // Loop over all valid combinations and build space points
+ for(Int_t icombo=0;icombo<ncombos;icombo++) {
+ THcDCHit* hits[4];
+ hits[0]=combos[icombo].pair1->hit1;
+ hits[1]=combos[icombo].pair1->hit2;
+ hits[2]=combos[icombo].pair2->hit1;
+ hits[3]=combos[icombo].pair2->hit2;
+ // Get Average Space point xt, yt
+ Double_t xt = combos[icombo].pair1->x + combos[icombo].pair2->x;
+ Double_t yt = combos[icombo].pair1->y + combos[icombo].pair2->y;
+ // Loop over space points
+ if(fNSpacePoints > 0) {
+ Int_t add_flag=1;
+ for(Int_t ispace=0;ispace<fNSpacePoints;ispace++) {
+ if(fSpacePoints[ispace].nhits > 0) {
+ Double_t sqdist_test = pow(xt - fSpacePoints[ispace].x,2) +
+ pow(yt - fSpacePoints[ispace].y,2);
+ // I (who is I) want to be careful if sqdist_test is bvetween 1 and
+ // 3 fSpacePointCriterion2. Let me ignore not add a new point the
+ if(sqdist_test < 3*fSpacePointCriterion2) {
+ add_flag = 0; // do not add a new space point
+ }
+ if(sqdist_test < fSpacePointCriterion2) {
+ // This is a real match
+ // Add the new hits to the existing space point
+ Int_t iflag[4];
+ iflag[0]=0;iflag[1]=0;iflag[2]=0;iflag[3]=0;
+ // Find out which of the four hits in the combo are already
+ // in the space point under consideration so that we don't
+ // add duplicate hits to the space point
+ for(Int_t isp_hit=0;isp_hit<fSpacePoints[ispace].nhits;isp_hit++) {
+ for(Int_t icm_hit=0;icm_hit<4;icm_hit++) { // Loop over combo hits
+ if(fSpacePoints[ispace].hits[isp_hit]==hits[icm_hit]) {
+ iflag[icm_hit] = 1;
+ }
+ }
+ }
+ // Remove duplicated pionts in the combo so we don't add
+ // duplicate hits to the space point
+ for(Int_t icm1=0;icm1<3;icm1++) {
+ for(Int_t icm2=icm1+1;icm2<4;icm2++) {
+ if(hits[icm1]==hits[icm2]) {
+ iflag[icm2] = 1;
+ }
+ }
+ }
+ // Add the unique combo hits to the space point
+ for(Int_t icm=0;icm<4;icm++) {
+ if(iflag[icm]==0) {
+ fSpacePoints[ispace].hits[fSpacePoints[ispace].nhits++] = hits[icm];
+ }
+ fSpacePoints[ispace].ncombos++;
+ }
+ }
+ }
+ }// End of loop over existing space points
+ // Create a new space point if more than 2*space_point_criteria
+ if(fNSpacePoints < MAX_SPACE_POINTS) {
+ if(add_flag) {
+ fSpacePoints[fNSpacePoints].nhits=2;
+ fSpacePoints[fNSpacePoints].ncombos=1;
+ fSpacePoints[fNSpacePoints].hits[0]=hits[0];
+ fSpacePoints[fNSpacePoints].hits[1]=hits[1];
+ fSpacePoints[fNSpacePoints].x = xt;
+ fSpacePoints[fNSpacePoints].y = yt;
+ if(hits[0] != hits[2] && hits[1] != hits[2]) {
+ fSpacePoints[fNSpacePoints].hits[fSpacePoints[fNSpacePoints].nhits++] = hits[2];
+ }
+ if(hits[0] != hits[3] && hits[1] != hits[3]) {
+ fSpacePoints[fNSpacePoints].hits[fSpacePoints[fNSpacePoints].nhits++] = hits[3];
+ }
+ fNSpacePoints++;
+ }
+ }
+ } else {// Create first space point
+ // This duplicates code above. Need to see if we can restructure
+ // to avoid
+ fSpacePoints[fNSpacePoints].nhits=2;
+ fSpacePoints[fNSpacePoints].ncombos=1;
+ fSpacePoints[fNSpacePoints].hits[0]=hits[0];
+ fSpacePoints[fNSpacePoints].hits[1]=hits[1];
+ fSpacePoints[fNSpacePoints].x = xt;
+ fSpacePoints[fNSpacePoints].y = yt;
+ if(hits[0] != hits[2] && hits[1] != hits[2]) {
+ fSpacePoints[fNSpacePoints].hits[fSpacePoints[fNSpacePoints].nhits++] = hits[2];
+ }
+ if(hits[0] != hits[3] && hits[1] != hits[3]) {
+ fSpacePoints[fNSpacePoints].hits[fSpacePoints[fNSpacePoints].nhits++] = hits[3];
+ }
+ fNSpacePoints++;
+ }//End check on 0 space points
+ }//End loop over combos
+ return(fNSpacePoints);
}
//_____________________________________________________________________________
-Int_t THcDriftChamber::Decode( const THaEvData& evdata )
+// HMS Specific?
+Int_t THcDriftChamber::DestroyPoorSpacePoints()
{
+ Int_t nhitsperplane[fNPlanes];
- // Get the Hall C style hitlist (fRawHitList) for this event
- Int_t nhits = THcHitList::DecodeToHitList(evdata);
+ Int_t spacepointsgood[fNSpacePoints];
+ Int_t ngood=0;
- // Let each plane get its hits
- Int_t nexthit = 0;
- for(Int_t ip=0;ip<fNPlanes;ip++) {
- nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
+ for(Int_t i=0;i<fNSpacePoints;i++) {
+ spacepointsgood[i] = 0;
+ }
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ Int_t nplanes_hit = 0;
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ nhitsperplane[ip] = 0;
+ }
+ // Count # hits in each plane for this space point
+ for(Int_t ihit=0;ihit<fSpacePoints[isp].nhits;ihit++) {
+ THcDCHit* hit=fSpacePoints[isp].hits[ihit];
+ // hit_order(hit) = ihit;
+ Int_t ip = hit->GetPlaneIndex();
+ nhitsperplane[ip]++;
+ }
+ // Count # planes that have hits
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if(nhitsperplane[ip] > 0) {
+ nplanes_hit++;
+ }
+ }
+ if(nplanes_hit >= fMinHits && nhitsperplane[YPlaneInd]>0
+ && nhitsperplane[YPlanePInd] > 0) {
+ spacepointsgood[ngood++] = isp; // Build list of good points
+ } else {
+ // cout << "Missing Y-hit!!";
+ }
}
-#if 0
- // fRawHitList is TClones array of THcRawDCHit objects
- for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
- THcRawDCHit* hit = (THcRawDCHit *) fRawHitList->At(ihit);
- // cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
- // << endl;
- for(Int_t imhit = 0; imhit < hit->fNHits; imhit++) {
- // cout << " " << imhit << " " << hit->fTDC[imhit]
- // << endl;
+ // Remove the bad space points
+ Int_t nremoved=fNSpacePoints-ngood;
+ fNSpacePoints = ngood;
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) { // New index num ber
+ Int_t osp=spacepointsgood[isp]; // Original index number
+ if(osp > isp) {
+ // Does this work, or do we have to copy each member?
+ // If it doesn't we should overload the = operator
+ fSpacePoints[isp] = fSpacePoints[osp];
}
}
- // cout << endl;
+ return nremoved;
+}
+
+//_____________________________________________________________________________
+// HMS Specific?
+ /*
+ Purpose and Methods : This routine loops over space points and
+ looks at all hits in the space
+ point. If more than 1 hit is in the same
+ plane then the space point is cloned with
+ all combinations of 1 wire per plane. The
+ requirements for cloning are: 1) at least
+ 4 planes fire, and 2) no more than 6 planes
+ have multiple hits.
+ */
+Int_t THcDriftChamber::SpacePointMultiWire()
+{
+ Int_t nhitsperplane[fNPlanes];
+ THcDCHit* hits_plane[fNPlanes][MAX_HITS_PER_POINT];
+
+ Int_t nsp_check;
+ Int_t nplanes_single;
+
+ Int_t nsp_tot=fNSpacePoints;
+
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ Int_t nplanes_hit = 0; // Number of planes with hits
+ Int_t nplanes_mult = 0; // Number of planes with multiple hits
+ Int_t nsp_new = 1;
+ Int_t newsp_num=0;
+
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ nhitsperplane[ip] = 0;
+ }
+ // Sort Space Points hits by plane
+ for(Int_t ihit=0;ihit<fSpacePoints[isp].nhits;ihit++) { // All hits in SP
+ THcDCHit* hit=fSpacePoints[isp].hits[ihit];
+ // hit_order Make a hash
+ // hash(hit) = ihit;
+ Int_t ip = hit->GetPlaneIndex();
+ hits_plane[ip][nhitsperplane[ip]++] = hit;
+ }
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ if(nhitsperplane[ip] > 0) {
+ nplanes_hit++;
+ nsp_new *= nhitsperplane[ip];
+ if(nhitsperplane[ip] > 1) nplanes_mult++;
+ }
+ }
+ --nsp_new;
+ nsp_check=nsp_tot + nsp_new;
+ nplanes_single = nplanes_hit - nplanes_mult;
+
+ // Check if cloning conditions are met
+ Int_t ntot = 0;
+ if(nplanes_hit >= 4 && nplanes_mult < 4 && nplanes_mult >0
+ && nsp_check < 20) {
+
+ // Order planes by decreasing # of hits
+
+ Int_t maxplane[fNPlanes];
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ maxplane[ip] = ip;
+ }
+ // Sort by decreasing # of hits
+ for(Int_t ip1=0;ip1<fNPlanes-1;ip1++) {
+ for(Int_t ip2=ip1+1;ip2<fNPlanes;ip2++) {
+ if(nhitsperplane[maxplane[ip2]] > nhitsperplane[maxplane[ip1]]) {
+ Int_t temp = maxplane[ip1];
+ maxplane[ip1] = maxplane[ip2];
+ maxplane[ip2] = temp;
+ }
+ }
+ }
+
+ // First fill clones with 1 hit each from the 3 planes with the most hits
+ for(Int_t n1=0;n1<nhitsperplane[maxplane[0]];n1++) {
+ for(Int_t n2=0;n2<nhitsperplane[maxplane[1]];n2++) {
+ for(Int_t n3=0;n3<nhitsperplane[maxplane[2]];n3++) {
+ ntot++;
+ newsp_num = nsp_tot + ntot - 2; // ntot will be 2 for first new
+ if(n1==0 && n2==0 && n3==0) newsp_num = isp; // Copy over original SP
+ fSpacePoints[newsp_num].x = fSpacePoints[isp].x;
+ fSpacePoints[newsp_num].y = fSpacePoints[isp].y;
+ fSpacePoints[newsp_num].nhits = nplanes_hit;
+ fSpacePoints[newsp_num].ncombos = fSpacePoints[isp].ncombos;
+ fSpacePoints[newsp_num].hits[0] = hits_plane[maxplane[0]][n1];
+ fSpacePoints[newsp_num].hits[1] = hits_plane[maxplane[1]][n2];
+ fSpacePoints[newsp_num].hits[2] = hits_plane[maxplane[2]][n3];
+ fSpacePoints[newsp_num].hits[3] = hits_plane[maxplane[3]][0];
+ if(nhitsperplane[maxplane[4]] == 1)
+ fSpacePoints[newsp_num].hits[4] = hits_plane[maxplane[4]][0];
+ if(nhitsperplane[maxplane[5]] == 1)
+ fSpacePoints[newsp_num].hits[5] = hits_plane[maxplane[5]][0];
+ }
+ }
+ }
+#if 0
+ // Loop over clones and order hits in the same way as parent SP
+ // Why do we have to order the hits.
+ for(Int_t i=0;i<ntot;i++) {
+ Int_t newsp_num= nsp_tot + i;
+ if(i == 1) newsp_num = isp;
+ for(Int_t j=0;j<nplanes_hit;j++) {
+ for(Int_t k=0;k<nplanes_hit;k++) {
+ THcDCHit* hit1 = fSpacePointHits[newsp_num][j];
+ THcDCHit* hit2 = fSpacePointHits[newsp_num][k];
+ if(hit_order(hit1) > hit_order(hit2)) {
+ THcDCHit* temp = fSpacePoints[newsp_num].hits[k];
+ fSpacePoints[newsp_num].hits[k] = fSpacePoints[newsp_num].hits[j];
+ fSpacePoints[newsp_num].hits[j] = temp;
+ }
+ }
+ }
+ }
#endif
+ nsp_tot += (ntot-1);
+ } else {
+ ntot=1;
+ }
+ }
+ assert (nsp_tot > 0);
+ Int_t nadded=0;
+ if(nsp_tot <= 20) {
+ nadded = nsp_tot - fNSpacePoints;
+ fNSpacePoints = nsp_tot;
+ }
- return nhits;
+ // In Fortran, fill in zeros.
+ return(nadded);
}
//_____________________________________________________________________________
-Int_t THcDriftChamber::ApplyCorrections( void )
+// HMS Specific?
+void THcDriftChamber::ChooseSingleHit()
{
- return(0);
+ // Look at all hits in a space point. If two hits are in the same plane,
+ // reject the one with the longer drift time.
+ Int_t goodhit[MAX_HITS_PER_POINT];
+
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ Int_t startnum = fSpacePoints[isp].nhits;
+
+ for(Int_t ihit=0;ihit<startnum;ihit++) {
+ goodhit[ihit] = 1;
+ }
+ // For each plane, mark all hits longer than the shortest drift time
+ for(Int_t ihit1=0;ihit1<startnum-1;ihit1++) {
+ THcDCHit* hit1 = fSpacePoints[isp].hits[ihit1];
+ Int_t plane1=hit1->GetPlaneIndex();
+ Double_t tdrift1 = hit1->GetTime();
+ for(Int_t ihit2=ihit1+1;ihit2<startnum;ihit2++) {
+ THcDCHit* hit2 = fSpacePoints[isp].hits[ihit2];
+ Int_t plane2=hit2->GetPlaneIndex();
+ Double_t tdrift2 = hit2->GetTime();
+ if(plane1 == plane2) {
+ if(tdrift1 > tdrift2) {
+ goodhit[ihit1] = 0;
+ } else {
+ goodhit[ihit2] = 0;
+ }
+ }
+ }
+ }
+ // Gather the remaining hits
+ Int_t finalnum = 0;
+ for(Int_t ihit=0;ihit<startnum;ihit++) {
+ if(goodhit[ihit] > 0) { // Keep this hit
+ if (ihit > finalnum) { // Move hit
+ fSpacePoints[isp].hits[finalnum++] = fSpacePoints[isp].hits[ihit];
+ } else {
+ finalnum++;
+ }
+ }
+ }
+ fSpacePoints[isp].nhits = finalnum;
+ }
}
+//_____________________________________________________________________________
+// Generic
+void THcDriftChamber::SelectSpacePoints()
+// This routine goes through the list of space_points and space_point_hits
+// found by find_space_points and only accepts those with
+// number of hits > min_hits
+// number of combinations > min_combos
+{
+ Int_t sp_count=0;
+ for(Int_t isp=0;isp<fNSpacePoints;isp++) {
+ // Include fEasySpacePoint because ncombos not filled in
+ if(fSpacePoints[isp].ncombos >= fMinCombos || fEasySpacePoint) {
+ if(fSpacePoints[isp].nhits >= fMinHits) {
+ fSpacePoints[sp_count++] = fSpacePoints[isp];
+ }
+ }
+ }
+ fNSpacePoints = sp_count;
+}
+
+/*
+*
+* Now we know rough hit positions in the chambers so we can make
+* wire velocity drift time corrections for each hit in the space point
+*
+* Assume all wires for a plane are read out on the same side (l/r or t/b).
+* If the wire is closer to horizontal, read out left/right. If nearer
+* vertical, assume top/bottom. (Note, this is not always true for the
+* SOS u and v planes. They have 1 card each on the side, but the overall
+* time offset per card will cancel much of the error caused by this. The
+* alternative is to check by card, rather than by plane and this is harder.
+*/
+
//_____________________________________________________________________________
-Int_t THcDriftChamber::CoarseTrack( TClonesArray& /* tracks */ )
+THcDriftChamber::~THcDriftChamber()
{
- // Calculation of coordinates of particle track cross point with scint
- // plane in the detector coordinate system. For this, parameters of track
- // reconstructed in THaVDC::CoarseTrack() are used.
- //
- // Apply corrections and reconstruct the complete hits.
- //
- // static const Double_t sqrt2 = TMath::Sqrt(2.);
-
- ApplyCorrections();
+ // Destructor. Remove variables from global list.
- return 0;
+ if( fIsSetup )
+ RemoveVariables();
+ if( fIsInit )
+ DeleteArrays();
+ if (fTrackProj) {
+ fTrackProj->Clear();
+ delete fTrackProj; fTrackProj = 0;
+ }
}
//_____________________________________________________________________________
-Int_t THcDriftChamber::FineTrack( TClonesArray& tracks )
+void THcDriftChamber::DeleteArrays()
{
- // Reconstruct coordinates of particle track cross point with scintillator
- // plane, and copy the data into the following local data structure:
- //
- // Units of measurements are meters.
+ // Delete member arrays. Used by destructor.
- // Calculation of coordinates of particle track cross point with scint
- // plane in the detector coordinate system. For this, parameters of track
- // reconstructed in THaVDC::FineTrack() are used.
+}
- return 0;
+//_____________________________________________________________________________
+inline
+void THcDriftChamber::Clear( const Option_t* )
+{
+ // Reset per-event data.
+ // fNhits = 0;
+
+ // fTrackProj->Clear();
+ fNhits = 0;
+
+}
+
+//_____________________________________________________________________________
+Int_t THcDriftChamber::ApplyCorrections( void )
+{
+ return(0);
}
ClassImp(THcDriftChamber)
diff --git a/src/THcDriftChamber.h b/src/THcDriftChamber.h
index a973c466a25804993d2e8c3576d05cefd9c1eb56..96c293c0fbe4c6a7600a712c3a907a804cc0413f 100644
--- a/src/THcDriftChamber.h
+++ b/src/THcDriftChamber.h
@@ -7,52 +7,40 @@
// //
///////////////////////////////////////////////////////////////////////////////
-#include "THaTrackingDetector.h"
-#include "THcHitList.h"
-#include "THcRawDCHit.h"
+#include "THaSubDetector.h"
#include "THcDriftChamberPlane.h"
-#include "TMath.h"
+#include "TClonesArray.h"
+
+#define MAX_SPACE_POINTS 50
+#define MAX_HITS_PER_POINT 20
+
+//#include "TMath.h"
//class THaScCalib;
class TClonesArray;
-class THcDriftChamber : public THaTrackingDetector, public THcHitList {
+class THcDriftChamber : public THaSubDetector {
public:
- THcDriftChamber( const char* name, const char* description = "",
- THaApparatus* a = NULL );
+ THcDriftChamber( const char* name, const char* description, Int_t chambernum,
+ THaDetectorBase* parent = NULL );
virtual ~THcDriftChamber();
- virtual Int_t Decode( const THaEvData& );
+ virtual Int_t Decode( const THaEvData& );
virtual EStatus Init( const TDatime& run_time );
- virtual Int_t CoarseTrack( TClonesArray& tracks );
- virtual Int_t FineTrack( TClonesArray& tracks );
+ virtual void AddPlane(THcDriftChamberPlane *plane);
virtual Int_t ApplyCorrections( void );
+ virtual void ProcessHits( void );
+ virtual Int_t FindSpacePoints( void ) ;
+
+ virtual void Clear( Option_t* opt="" );
// Int_t GetNHits() const { return fNhit; }
Int_t GetNTracks() const { return fTrackProj->GetLast()+1; }
const TClonesArray* GetTrackHits() const { return fTrackProj; }
- Int_t GetNWires(Int_t plane) const { return fNWires[plane-1];}
- Int_t GetNChamber(Int_t plane) const { return fNChamber[plane-1];}
- Int_t GetWireOrder(Int_t plane) const { return fWireOrder[plane-1];}
- Int_t GetPitch(Int_t plane) const { return fPitch[plane-1];}
- Int_t GetCentralWire(Int_t plane) const { return fCentralWire[plane-1];}
- Int_t GetTdcWinMin(Int_t plane) const { return fTdcWinMin[plane-1];}
- Int_t GetTdcWinMax(Int_t plane) const { return fTdcWinMax[plane-1];}
-
- Double_t GetPlaneTimeZero(Int_t plane) const { return fPlaneTimeZero[plane-1];}
-
- Double_t GetNSperChan() const { return fNSperChan;}
-
- Double_t GetCenter(Int_t plane) const {
- Int_t chamber = GetNChamber(plane)-1;
- return
- fXCenter[chamber]*sin(fAlphaAngle[plane-1]) +
- fYCenter[chamber]*cos(fAlphaAngle[plane-1]);
- }
// friend class THaScCalib;
THcDriftChamber(); // for ROOT I/O
@@ -61,65 +49,58 @@ protected:
// Calibration
// Per-event data
+ Int_t fNhits;
- // Potential Hall C parameters. Mostly here for demonstration
- Int_t fNPlanes;
- char** fPlaneNames;
- Int_t fNChambers;
-
- Double_t fNSperChan; /* TDC bin size */
- Double_t fWireVelocity;
-
- // Each of these will be dimensioned with the number of chambers
- Double_t* fXCenter;
- Double_t* fYCenter;
-
- // Each of these will be dimensioned with the number of planes
- // A THcDriftChamberPlane class object will need to access the value for
- // its plane number. Should we have a Get method for each or
- Int_t* fTdcWinMin;
- Int_t* fTdcWinMax;
- Int_t* fCentralTime;
- Int_t* fNWires; // Number of wires per plane
- Int_t* fNChamber;
- Int_t* fWireOrder;
- Int_t* fDriftTimeSign;
-
- Double_t* fZPos;
- Double_t* fAlphaAngle;
- Double_t* fBetaAngle;
- Double_t* fGammaAngle;
- Double_t* fPitch;
- Double_t* fCentralWire;
- Double_t* fPlaneTimeZero;
-
- THcDriftChamberPlane** fPlanes; // List of plane objects
+ Int_t fNPlanes; // Number of planes in the chamber
+
+ Int_t fChamberNum;
+
+ // HMS Specific
+ Int_t YPlaneInd; // Index of Yplane for this chamber
+ Int_t YPlanePInd; // Index of Yplanep for this chamber
+ Int_t YPlaneNum; // Absolute plane number of Yplane
+ Int_t YPlanePNum; // Absolute plane number of Yplanep
+
+ // Parameters
+ Int_t fMinHits; // Minimum hits required to do something
+ Int_t fMaxHits; // Maximum required to do something
+ Int_t fMinCombos; // Minimum # pairs in a space point
+ Int_t fRemove_Sppt_If_One_YPlane;
+
+ Double_t fXCenter;
+ Double_t fYCenter;
+ Double_t fSpacePointCriterion;
+ Double_t fSpacePointCriterion2;
+
+ THcDriftChamberPlane* fPlanes[20]; // List of plane 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 ClearEvent();
void DeleteArrays();
virtual Int_t ReadDatabase( const TDatime& date );
virtual Int_t DefineVariables( EMode mode = kDefine );
void Setup(const char* name, const char* description);
+ Int_t FindEasySpacePoint(Int_t yplane_hitind, Int_t yplanep_hitind);
+ Int_t FindHardSpacePoints(void);
+ Int_t DestroyPoorSpacePoints(void);
+ Int_t SpacePointMultiWire(void);
+ void ChooseSingleHit(void);
+ void SelectSpacePoints(void);
+
+ THcDCHit* fHits[10000]; /* All hits for this chamber */
+ // A simple structure until we figure out what we are doing.
+ struct SpacePoint {
+ Double_t x;
+ Double_t y;
+ Int_t nhits;
+ Int_t ncombos;
+ THcDCHit* hits[MAX_HITS_PER_POINT];
+ };
+ SpacePoint fSpacePoints[MAX_SPACE_POINTS];
+ Int_t fNSpacePoints;
+ Int_t fEasySpacePoint; /* This event is an easy space point */
ClassDef(THcDriftChamber,0) // Drift Chamber class
};
diff --git a/src/THcDriftChamberPlane.cxx b/src/THcDriftChamberPlane.cxx
index 816d27a650ea46865c1bd42f0655f5751de1ae1d..b23cfe21c160216f272ab7c3abd6cf807802dc2d 100644
--- a/src/THcDriftChamberPlane.cxx
+++ b/src/THcDriftChamberPlane.cxx
@@ -14,7 +14,7 @@
#include "THcGlobals.h"
#include "THcParmList.h"
#include "THcHitList.h"
-#include "THcDriftChamber.h"
+#include "THcDC.h"
#include "THcHodoscope.h"
#include "TClass.h"
@@ -97,11 +97,11 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
gHcParms->LoadParmValues((DBRequest*)&list,prefix);
// Retrieve parameters we need from parent class
- THcDriftChamber* fParent;
+ THcDC* fParent;
- fParent = (THcDriftChamber*) GetParent();
+ fParent = (THcDC*) GetParent();
// These are single variables here, but arrays in THcDriftChamber.
- fNChamber = fParent->GetNChamber(fPlaneNum);
+ fChamberNum = fParent->GetNChamber(fPlaneNum);
fNWires = fParent->GetNWires(fPlaneNum);
fWireOrder = fParent->GetWireOrder(fPlaneNum);
fPitch = fParent->GetPitch(fPlaneNum);
@@ -113,7 +113,34 @@ Int_t THcDriftChamberPlane::ReadDatabase( const TDatime& date )
fNSperChan = fParent->GetNSperChan();
- cout << fPlaneNum << " " << fNWires << endl;
+ // 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 alpha = fParent->GetAlphaAngle(fPlaneNum);
+ Double_t beta = fParent->GetBetaAngle(fPlaneNum);
+ Double_t gamma = fParent->GetGammaAngle(fPlaneNum);
+ Double_t cosalpha = TMath::Cos(alpha);
+ Double_t sinalpha = TMath::Sin(alpha);
+ Double_t cosbeta = TMath::Cos(beta);
+ 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 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;
+
+ cout << fPlaneNum << " " << fNWires << " " << fWireOrder << endl;
fTTDConv = new THcDCLookupTTDConv(DriftMapFirstBin,fPitch/2,DriftMapBinSize,
NumDriftMapBins,DriftMap);
@@ -242,7 +269,7 @@ Int_t THcDriftChamberPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
- rawtdc*fNSperChan + fPlaneTimeZero;
// How do we get this start time from the hodoscope to here
// (or at least have it ready by coarse process)
- new( (*fHits)[nextHit++] ) THcDCHit(wire, rawtdc, time);
+ new( (*fHits)[nextHit++] ) THcDCHit(wire, rawtdc, time, this);
}
wire_last = wireNum;
}
diff --git a/src/THcDriftChamberPlane.h b/src/THcDriftChamberPlane.h
index 7307a4acf649f75679223878d27d2607747c3798..78faf37fefb33635e0c8ac1546c2ea9caf34abe2 100644
--- a/src/THcDriftChamberPlane.h
+++ b/src/THcDriftChamberPlane.h
@@ -48,6 +48,16 @@ class THcDriftChamberPlane : public THaSubDetector {
{ assert( i>=1 && i<=GetNWires() );
return (THcDCWire*)fWires->UncheckedAt(i-1); }
+ Int_t GetNHits() const { return fHits->GetLast()+1; }
+ TClonesArray* GetHits() const { return fHits; }
+
+ Int_t GetPlaneNum() const { return fPlaneNum; }
+ Int_t GetChamberNum() const { return fChamberNum; }
+ void SetPlaneIndex(Int_t index) { fPlaneIndex = index; }
+ Int_t GetPlaneIndex() { return fPlaneIndex; }
+ Double_t GetXsp() const { return fXsp; }
+ Double_t GetYsp() const { return fYsp; }
+
protected:
TClonesArray* fParentHitList;
@@ -56,7 +66,8 @@ class THcDriftChamberPlane : public THaSubDetector {
TClonesArray* fWires;
Int_t fPlaneNum;
- Int_t fNChamber;
+ Int_t fPlaneIndex; /* Index of this plane within it's chamber */
+ Int_t fChamberNum;
Int_t fNWires;
Int_t fWireOrder;
Int_t fTdcWinMin;
@@ -64,6 +75,8 @@ class THcDriftChamberPlane : public THaSubDetector {
Double_t fPitch;
Double_t fCentralWire;
Double_t fPlaneTimeZero;
+ Double_t fXsp;
+ Double_t fYsp;
Double_t fCenter;