diff --git a/Makefile b/Makefile
index 6395c76eccf20788b53909404cd2801d02d16efb..83f6afeea8f44933f070a53c101b5673fcee8b85 100644
--- a/Makefile
+++ b/Makefile
@@ -28,7 +28,8 @@ SRC = src/THcInterface.cxx src/THcParmList.cxx src/THcAnalyzer.cxx \
src/THcRaster.cxx\
src/THcRasteredBeam.cxx\
src/THcRasterRawHit.cxx \
- src/THcScalerEvtHandler.cxx
+ src/THcScalerEvtHandler.cxx \
+ src/THcHodoEff.cxx
# Name of your package.
# The shared library that will be built will get the name lib$(PACKAGE).so
diff --git a/examples/hodtest.C b/examples/hodtest.C
index 9a09d917c57cc757ae4bf717e62b81826ac6d5bb..aec31be3f2cb021bf49708eac5c6d0326dae4818 100644
--- a/examples/hodtest.C
+++ b/examples/hodtest.C
@@ -36,7 +36,8 @@
gHaApps->Add( HMS );
// Add hodoscope
- HMS->AddDetector( new THcHodoscope("hod", "Hodoscope" ));
+ THcHodoscope* hms_hodoscope = new THcHodoscope("hod","Hodoscope");
+ HMS->AddDetector( hms_hodoscope );
HMS->AddDetector( new THcShower("cal", "Shower" ));
HMS->AddDetector( new THcDC("dc", "Drift Chambers" ));
THcAerogel* aerogel = new THcAerogel("aero", "Aerogel Cerenkov" );
@@ -47,10 +48,15 @@
THaApparatus* SOS = new THcHallCSpectrometer("S","SOS");
gHaApps->Add( SOS );
// Add detectors
- SOS->AddDetector( new THcHodoscope("hod", "Hodoscope" ));
+ THcHodoscope* sos_hodoscope = new THcHodoscope("hod","Hodoscope");
+ SOS->AddDetector( sos_hodoscope);
SOS->AddDetector( new THcShower("cal", "Shower" ));
SOS->AddDetector( new THcDC("dc", "Drift Chambers" ));
+ gHaPhysics->Add(new THcHodoEff("hhodeff","HMS Hodoscope Efficiencies","H.hod"));
+ gHaPhysics->Add(new THcHodoEff("shodeff","SOS Hodoscope Efficiencies","S.hod"));
+
+
// Set up the analyzer - we use the standard one,
// but this could be an experiment-specific one as well.
// The Analyzer controls the reading of the data, executes
diff --git a/src/HallC_LinkDef.h b/src/HallC_LinkDef.h
index 39974874083e9012ee8a36cb0b30853f5a76a8e5..93be610ca10cea679e3a37d2e11332214c8782f3 100644
--- a/src/HallC_LinkDef.h
+++ b/src/HallC_LinkDef.h
@@ -53,5 +53,6 @@
#pragma link C++ class THcRasteredBeam+;
#pragma link C++ class THcRasterRawHit+;
#pragma link C++ class THcScalerEvtHandler+;
+#pragma link C++ class THcHodoEff+;
#endif
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index d714583e0b63a73e83acdb681426e45a19b4a42e..7c1216dd2373a90f3d69b613dd1f07ce4de6c3e8 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -529,6 +529,8 @@ Int_t THcDC::CoarseTrack( TClonesArray& tracks )
// CalcFocalPlaneCoords. Aren't our tracks already in focal plane coords
// We should have some kind of track ID so that the THaTrack can be
// associate back with the DC track
+ // Assign the track number
+ theTrack->SetTrkNum(itrack+1);
}
}
diff --git a/src/THcHodoEff.cxx b/src/THcHodoEff.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..bff2a480207a42039fc2e6c2f5d9f93fbb5815ed
--- /dev/null
+++ b/src/THcHodoEff.cxx
@@ -0,0 +1,354 @@
+///////////////////////////////////////////////////////////////////////////////
+// //
+// THcHodoEff //
+// //
+// Class for accumulating statistics for and calculating hodoscope //
+// efficiencies. //
+// //
+// Moddled after VDCeff //
+// For now trying to emulate work done in h_scin_eff/h_scin_eff_shutdown //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THcHodoEff.h"
+#include "THaApparatus.h"
+#include "THcHodoHit.h"
+#include "THcGlobals.h"
+#include "THcParmList.h"
+
+using namespace std;
+
+//_____________________________________________________________________________
+THcHodoEff::THcHodoEff (const char *name, const char* description,
+ const char* hodname) :
+ THaPhysicsModule(name, description), fName(hodname), fHod(NULL), fNevt(0)
+{
+
+}
+
+//_____________________________________________________________________________
+THcHodoEff::~THcHodoEff()
+{
+ // Destructor
+
+ RemoveVariables();
+}
+
+//_____________________________________________________________________________
+void THcHodoEff::Reset( Option_t* opt )
+// Clear event-by-event data
+{
+ Clear(opt);
+}
+
+//_____________________________________________________________________________
+Int_t THcHodoEff::Begin( THaRunBase* )
+{
+ // Start of analysis
+
+ if (!IsOK() ) return -1;
+
+ // Book any special histograms here
+
+ fNevt = 0;
+
+ // Clear all the accumulators here
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ fHitPlane[ip] = 0;
+ for(Int_t ic=0;ic<fNCounters[ip];ic++) {
+ fStatPosHit[ip][ic] = 0;
+ fStatNegHit[ip][ic] = 0;
+ fStatAndHit[ip][ic] = 0;
+ fStatOrHit[ip][ic] = 0;
+ for(Int_t idel=0;idel<20;idel++) {
+ fStatTrkDel[ip][ic][idel] = 0;
+ fStatAndHitDel[ip][ic][idel] = 0;
+ }
+ }
+ }
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcHodoEff::End( THaRunBase* )
+{
+ // End of analysis
+
+ return 0;
+}
+
+
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcHodoEff::Init( const TDatime& run_time )
+{
+ // Initialize THcHodoEff physics module
+
+ // const char* const here = "Init";
+
+ // Standard initialization. Calls ReadDatabase(), ReadRunDatabase(),
+ // and DefineVariables() (see THaAnalysisObject::Init)
+
+ fHod = dynamic_cast<THcHodoscope*>
+ ( FindModule( fName.Data(), "THcHodoscope"));
+
+ fSpectro = static_cast<THaSpectrometer*>(fHod->GetApparatus());
+
+ if( THaPhysicsModule::Init( run_time ) != kOK )
+ return fStatus;
+
+ cout << "THcHodoEff::Init nplanes=" << fHod->GetNPlanes() << endl;
+ cout << "THcHodoEff::Init Apparatus = " << fHod->GetName() <<
+ " " <<
+ (fHod->GetApparatus())->GetName() << endl;
+
+ return fStatus = kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcHodoEff::Process( const THaEvData& evdata )
+{
+ // Accumulate statistics for efficiency
+
+ // const char* const here = "Process";
+
+ if( !IsOK() ) return -1;
+
+ // Project the golden track to each
+ // plane. Need to get track at Focal Plane, not tgt.
+ //
+ // Assumes that planes are X, Y, X, Y
+ THaTrack* theTrack = fSpectro->GetGoldenTrack();
+ Int_t trackIndex = theTrack->GetTrkNum()-1;
+ // Since fSpectro knows the index of the golden track, we can
+ // get other information about the track from fSpectro.
+ // Need to remove the specialized stuff from fGoldenTrack
+
+ if(!theTrack) return 0;
+
+ // May make these member variables
+ Double_t hitPos[fNPlanes];
+ Double_t hitDistance[fNPlanes];
+ Int_t hitCounter[fNPlanes];
+ Int_t checkHit[fNPlanes];
+ Bool_t goodTdcBothSides[fNPlanes];
+ Bool_t goodTdcOneSide[fNPlanes];
+
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ // Should really have plane object self identify as X or Y
+ if(ip%2 == 0) { // X Plane
+ hitPos[ip] = theTrack->GetX() + theTrack->GetTheta()*fPosZ[ip];
+ } else { // Y Plane
+ hitPos[ip] = theTrack->GetY() + theTrack->GetPhi()*fPosZ[ip];
+ }
+ // Which counter does track pass through?
+ hitCounter[ip] = TMath::Min(TMath::Max(TMath::Nint((hitPos[ip] - fCenterFirst[ip])/fSpacing[ip]) + 1,1),fNCounters[ip]);
+ // How far from paddle center is track?
+ if(ip%2 == 0) { // X Plane
+ hitDistance[ip] = hitPos[ip] - (fSpacing[ip]*(hitCounter[ip]-1) +
+ fCenterFirst[ip]);
+ } else { // Y Plane
+ hitDistance[ip] = hitPos[ip] - (fSpacing[ip]*(hitCounter[ip]-1) -
+ fCenterFirst[ip]);
+ }
+ }
+
+ // Fill dpos histograms and set checkHit for each plane.
+ // dpos stuff not implemented
+ // Why do dpos stuff here, does any other part need the dpos historgrams
+ // Look to VDCEff code to see how to create and fill histograms
+
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ Int_t hitcounter=hitCounter[ip];
+ goodTdcBothSides[ip] = kFALSE;
+ goodTdcOneSide[ip] = kFALSE;
+ checkHit[ip] = 2;
+ Int_t nphits=fPlanes[ip]->GetNScinHits();
+ TClonesArray* hodoHits = fPlanes[ip]->GetHits();
+ for(Int_t ihit=0;ihit<nphits;ihit++) {
+ THcHodoHit* hit=(THcHodoHit*) hodoHits->At(ihit);
+ Int_t counter = hit->GetPaddleNumber();
+ if(counter == hitcounter) {
+ checkHit[ip] = 0;
+ } else {
+ if(TMath::Abs(counter-hitcounter) == 1 && checkHit[ip] != 0) {
+ checkHit[ip] = 1;
+ }
+ }
+ }
+ }
+
+ // Record position differences between track and center of scin
+ // and increment 'should have hit' counters
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ Int_t hitcounter = hitCounter[ip];
+ Double_t dist = hitDistance[ip];
+ if(TMath::Abs(dist) <= fStatSlop &&
+ theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq) {
+ // && hsshtrk >= 0.05
+ Double_t delta = theTrack->GetDp();
+ Int_t idel = TMath::Floor(delta+10.0);
+ // Should
+ if(idel >=0 && idel < 20) {
+ fStatTrkDel[ip][hitcounter][idel]++;
+ }
+ // lookat[ip] = TRUE;
+ }
+ fHitPlane[ip] = 0;
+ }
+ // Is there a hit on or adjacent to paddle that track
+ // passes through?
+
+ // May collapse this loop into last
+
+ // record the hits as a "didhit" if track is near center of
+ // scintillator, the chisqared of the track is good and it is the
+ // first "didhit" in that plane.
+
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ Int_t hitcounter = hitCounter[ip];
+ Double_t dist = hitDistance[ip];
+ Int_t nphits=fPlanes[ip]->GetNScinHits();
+ TClonesArray* hodoHits = fPlanes[ip]->GetHits();
+ for(Int_t ihit=0;ihit<nphits;ihit++) {
+ THcHodoHit* hit=(THcHodoHit*) hodoHits->At(ihit);
+ Int_t counter = hit->GetPaddleNumber();
+ // Finds first best hit
+ Bool_t onTrack, goodScinTime, goodTdcNeg, goodTdcPos;
+ fHod->GetFlags(trackIndex,ip,ihit,
+ onTrack, goodScinTime, goodTdcNeg, goodTdcPos);
+ if(TMath::Abs(dist) <= fStatSlop &&
+ TMath::Abs(hitcounter-counter) <= checkHit[ip] &&
+ fHitPlane[ip] == 0 &&
+ theTrack->GetChi2()/theTrack->GetNDoF() <= fMaxChisq) {
+ // && hsshtrk >= 0.05
+ fHitPlane[ip]++;
+
+ // Need to find out hgood_tdc_pos(igoldentrack,ihit) and neg
+ if(goodTdcPos) {
+ if(goodTdcNeg) { // Both fired
+ fStatPosHit[ip][hitcounter]++;
+ fStatNegHit[ip][hitcounter]++;
+ fStatAndHit[ip][hitcounter]++;
+ fStatOrHit[ip][hitcounter]++;
+ Double_t delta = theTrack->GetDp();
+ Int_t idel = TMath::Floor(delta+10.0);
+ if(idel >=0 && idel < 20) {
+ fStatAndHitDel[ip][hitcounter][idel]++;
+ }
+ } else {
+ fStatPosHit[ip][hitcounter]++;
+ fStatOrHit[ip][hitcounter]++;
+ }
+ } else if(goodTdcNeg) {
+ fStatNegHit[ip][hitcounter]++;
+ fStatOrHit[ip][hitcounter]++;
+ }
+ }
+ // Increment pos/neg/both fired. Track independent, so
+ // no chisquared cut, but note that only scintillators on the
+ // track are examined.
+ if(goodTdcPos) {
+ if(goodTdcNeg) {
+ fBothGood[ip][hitcounter]++;
+ } else {
+ fPosGood[ip][hitcounter]++;
+ }
+ } else if (goodTdcNeg) {
+ fNegGood[ip][hitcounter]++;
+ }
+ // Determine if one or both PMTs had a good tdc
+
+ if(goodTdcPos && goodTdcNeg) {
+ goodTdcBothSides[ip] = kTRUE;
+ }
+ if(goodTdcPos || goodTdcNeg) {
+ goodTdcOneSide[ip] = kTRUE;
+ }
+ }
+ }
+ /*
+ For each plane, see of other 3 fired. This means that they were enough
+ to form a 3/4 trigger, and so the fraction of times this plane fired is
+ the plane trigger efficiency. NOTE: we only require a TDC hit, not a
+ TDC hit within the SCIN 3/4 trigger window, so high rates will make
+ this seem better than it is. Also, make sure we're not near the edge
+ of the hodoscope (at the last plane), using the same hhodo_slop param. as for h_tof.f
+ NOTE ALSO: to make this check simpler, we are assuming that all planes
+ have identical active areas. y_scin = y_cent + y_offset, so shift track
+ position by offset for comparing to edges.
+ */
+
+ // Need to add calculation and cuts on
+ // xatback and yatback in order to set the
+ // htrig_hododidflag, htrig_hodoshouldflag and otherthreehit flags
+ //
+
+ ++fNevt;
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcHodoEff::ReadDatabase( const TDatime& date )
+{
+ // Read database. Gets variable needed for efficiency calculation
+ // Get # of planes and their z positions here.
+
+ fNPlanes = fHod->GetNPlanes();
+ fPlanes = new THcScintillatorPlane* [fNPlanes];
+ fPosZ = new Double_t[fNPlanes];
+ fSpacing = new Double_t[fNPlanes];
+ fCenterFirst = new Double_t[fNPlanes];
+ fNCounters = new Int_t[fNPlanes];
+ fHodoSlop = new Double_t[fNPlanes];
+
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ fPlanes[ip] = fHod->GetPlane(ip);
+ fPosZ[ip] = fPlanes[ip]->GetZpos() + 0.5*fPlanes[ip]->GetDzpos();
+ fSpacing[ip] = fPlanes[ip]->GetSpacing();\
+ fCenterFirst[ip] = fPlanes[ip]->GetPosCenter(0);
+ fNCounters[ip] = fPlanes[ip]->GetNelem();
+ }
+
+ char prefix[2];
+ prefix[0] = tolower((fHod->GetApparatus())->GetName()[0]);
+ prefix[1] = '\0';
+
+ DBRequest list[]={
+ {"stat_slop", &fStatSlop, kDouble},
+ {"stat_maxchisq",&fMaxChisq, kDouble},
+ {"hodo_slop", fHodoSlop, kDouble, fNPlanes},
+ {0}
+ };
+ // fMaxShTrk = 0.05; // For cut on fraction of momentum seen in shower
+ gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+
+ cout << "THcHodoEff::ReadDatabase nplanes=" << fHod->GetNPlanes() << endl;
+ // Setup statistics arrays
+ // Better method to put this in?
+ // These all need to be cleared in Begin
+ fHitPlane = new Int_t[fNPlanes];
+ fStatTrkDel.resize(fNPlanes);
+ fStatAndHitDel.resize(fNPlanes);
+ fStatPosHit.resize(fNPlanes);
+ fStatNegHit.resize(fNPlanes);
+ fStatAndHit.resize(fNPlanes);
+ fStatOrHit.resize(fNPlanes);
+ for(Int_t ip=0;ip<fNPlanes;ip++) {
+ fStatTrkDel[ip].resize(fNCounters[ip]);
+ fStatAndHitDel[ip].resize(fNCounters[ip]);
+ fStatPosHit[ip].resize(fNCounters[ip]);
+ fStatNegHit[ip].resize(fNCounters[ip]);
+ fStatAndHit[ip].resize(fNCounters[ip]);
+ fStatOrHit[ip].resize(fNCounters[ip]);
+ for(Int_t ic=0;ic<fNCounters[ip];ic++) {
+ fStatTrkDel[ip][ic].resize(20); // Max this settable
+ fStatAndHitDel[ip][ic].resize(20); // Max this settable
+ }
+ }
+ return kOK;
+}
+//_____________________________________________________________________________
+ClassImp(THcHodoEff)
+////////////////////////////////////////////////////////////////////////////////
+
diff --git a/src/THcHodoEff.h b/src/THcHodoEff.h
new file mode 100644
index 0000000000000000000000000000000000000000..eaed7f07f95f7e23365abfbde41428daffd59251
--- /dev/null
+++ b/src/THcHodoEff.h
@@ -0,0 +1,72 @@
+#ifndef ROOT_THcHodoEff
+#define ROOT_THcHodoEff
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// THcHodoEff //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THaPhysicsModule.h"
+#include "THcHodoscope.h"
+#include "THaSpectrometer.h"
+#include "THaTrack.h"
+
+class THcHodoEff : public THaPhysicsModule {
+public:
+ THcHodoEff( const char* name, const char* description, const char* hodname);
+ virtual ~THcHodoEff();
+
+ virtual Int_t Begin( THaRunBase* r=0 );
+ virtual Int_t End( THaRunBase* r=0 );
+ virtual EStatus Init( const TDatime& run_time );
+ virtual Int_t Process( const THaEvData& );
+
+ void Reset( Option_t* opt="" );
+
+protected:
+
+ virtual Int_t ReadDatabase( const TDatime& date);
+
+ // Data needed for efficiency calculation for one Hodoscope paddle
+
+ Double_t* fZPos; //
+
+ TString fName; // Name of hodoscope
+ THcHodoscope* fHod; // Hodscope object
+ THaSpectrometer* fSpectro; // Spectrometer object
+
+ Long64_t fNevt;
+
+ // Information about the hodoscopes that we get from the
+ // THcHodoscope object
+
+ Int_t fNPlanes;
+ THcScintillatorPlane** fPlanes;
+ Double_t* fPosZ;
+ Double_t* fSpacing;
+ Double_t* fCenterFirst;
+ Int_t* fNCounters;
+ // Int_t fMaxNcounters;
+ Double_t fStatSlop;
+ Double_t fMaxChisq;
+ Double_t* fHodoSlop;
+
+ // Arrays for accumulating statistics
+ vector<vector<vector<Int_t> > > fHitShould;
+ vector<vector<vector<Int_t> > > fStatAndHitDel;
+ vector<vector<vector<Int_t> > > fStatTrkDel;
+ vector<vector<Int_t> > fStatPosHit;
+ vector<vector<Int_t> > fStatNegHit;
+ vector<vector<Int_t> > fStatAndHit;
+ vector<vector<Int_t> > fStatOrHit;
+ vector<vector<Int_t> > fBothGood;
+ vector<vector<Int_t> > fNegGood;
+ vector<vector<Int_t> > fPosGood;
+
+ Int_t* fHitPlane;
+
+ ClassDef(THcHodoEff,0) // Hodoscope efficiency module
+};
+
+#endif