diff --git a/src/Scandalizer.cxx b/src/Scandalizer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9d802a413a9bd6cf25a491603150ec0cb373f4df
--- /dev/null
+++ b/src/Scandalizer.cxx
@@ -0,0 +1,312 @@
+#include "Scandalizer.h"
+
+using namespace std;
+
+namespace hcana {
+
+Int_t Scandalizer::ReadOneEvent()
+{
+ // Read one event from current run (fRun) and raw-decode it using the
+ // current decoder (fEvData)
+
+ if( fDoBench ) fBench->Begin("RawDecode");
+
+ bool to_read_file = false;
+ if( !fEvData->IsMultiBlockMode() ||
+ (fEvData->IsMultiBlockMode() && fEvData->BlockIsDone()) )
+ to_read_file = true;
+
+ // Find next event buffer in CODA file. Quit if error.
+ Int_t status = THaRunBase::READ_OK;
+ if (to_read_file)
+ status = fRun->ReadEvent();
+
+ // there may be a better place to do this, but this works
+ if (fWantCodaVers > 0) {
+ fEvData->SetCodaVersion(fWantCodaVers);
+ } else {
+ fEvData->SetCodaVersion(fRun->GetCodaVersion());
+ }
+
+ switch( status ) {
+ case THaRunBase::READ_OK:
+ // Decode the event
+ if (to_read_file) {
+ status = fEvData->LoadEvent( fRun->GetEvBuffer() );
+ } else {
+ status = fEvData->LoadFromMultiBlock( ); // load next event in block
+ }
+ switch( status ) {
+ case THaEvData::HED_WARN:
+ //std::cout << "HED_WARN\n";
+ case THaEvData::HED_OK: // fall through
+ status = THaRunBase::READ_OK;
+ Incr(kNevRead);
+ break;
+ case THaEvData::HED_ERR:
+ // Decoding error
+ status = THaRunBase::READ_ERROR;
+ Incr(kDecodeErr);
+ break;
+ case THaEvData::HED_FATAL:
+ status = THaRunBase::READ_FATAL;
+ break;
+ }
+ break;
+
+ case THaRunBase::READ_EOF: // fall through
+ std::cout << " READ_EOF\n";
+ case THaRunBase::READ_FATAL:
+ // Just exit on EOF - don't count it
+ break;
+ default:
+ Incr(kCodaErr);
+ break;
+ }
+
+ if( fDoBench ) fBench->Stop("RawDecode");
+ return status;
+}
+
+Int_t Scandalizer::Process( THaRunBase* run )
+{
+ // Process the given run. Loop over all events in the event range and
+ // analyze all apparatuses defined in the global apparatus list.
+ // Fill Event structure if it is defined.
+ // If Event and Filename are defined, then fill the output tree with Event
+ // and write the file.
+
+ static const char* const here = "Process";
+
+ if( !run ) {
+ if( fRun )
+ run = fRun;
+ else
+ return -1;
+ }
+
+ //--- Initialization. Creates fFile, fOutput, and fEvent if necessary.
+ // Also copies run to fRun if run is different from fRun
+ Int_t status = Init( run );
+ if( status != 0 ) {
+ return status;
+ }
+
+ // Restart "Total" since it is stopped in Init()
+ fBench->Begin("Total");
+
+ //--- Re-open the data source. Should succeed since this was tested in Init().
+ if( (status = fRun->Open()) != THaRunBase::READ_OK ) {
+ Error( here, "Failed to re-open the input file. "
+ "Make sure the file still exists.");
+ fBench->Stop("Total");
+ return -4;
+ }
+
+ // Make the current run available globally - the run parameters are
+ // needed by some modules
+ gHaRun = fRun;
+
+ // Enable/disable helicity decoding as requested
+ fEvData->EnableHelicity( HelicityEnabled() );
+ // Set decoder reporting level. FIXME: update when THaEvData is updated
+ fEvData->SetVerbose( (fVerbose>2) );
+ fEvData->SetDebug( (fVerbose>3) );
+
+ // Informational messages
+ if( fVerbose>1 ) {
+ cout << "Decoder: helicity "
+ << (fEvData->HelicityEnabled() ? "enabled" : "disabled")
+ << endl;
+ cout << endl << "Starting analysis" << endl;
+ }
+ if( fVerbose>2 && fRun->GetFirstEvent()>1 )
+ cout << "Skipping " << fRun->GetFirstEvent() << " events" << endl;
+
+ //--- The main event loop.
+
+ fNev = 0;
+ bool terminate = false, fatal = false;
+ UInt_t nlast = fRun->GetLastEvent();
+ fAnalysisStarted = kTRUE;
+ BeginAnalysis();
+ if( fFile ) {
+ fFile->cd();
+ fRun->Write("Run_Data"); // Save run data to first ROOT file
+ }
+ void (*prev_handler)(int);
+ prev_handler = signal (SIGINT, handle_sig);
+
+ while ( !terminate && (fNev < nlast)) {
+ //&& (status = ReadOneEvent()) != THaRunBase::READ_EOF )
+ //std::cout << " evtype(last) " << fEvData->GetEvType() << "\n";
+ status = ReadOneEvent();
+ //std::cout << " status " << status << "\n";
+ //std::cout << " evtype " << fEvData->GetEvType() << "\n";
+
+ // If an interupt signal is sent (ctrl-c)
+ // //
+ if(sig_caught) {
+ terminate = true;
+ break;
+ }
+
+
+ if( status == THaRunBase::READ_EOF ) {
+ //if( fEvData->GetEvType() != 20 ) {
+ // std::cout << " sleeping....";
+ // //gSystem->Sleep(1000);
+ //}else {
+ //}
+ break;
+ }
+
+ //--- Skip events with errors, unless fatal
+ if( status == THaRunBase::READ_FATAL ){
+ std::cout << " READ_FATAL\n";
+ break;
+ }
+ if( status != THaRunBase::READ_OK ){
+ continue;
+ }
+
+ UInt_t evnum = fEvData->GetEvNum();
+
+ // Count events according to the requested mode
+ // Whether or not to ignore events prior to fRun->GetFirstEvent()
+ // is up to the analysis routines.
+ switch(fCountMode) {
+ case kCountPhysics:
+ if( fEvData->IsPhysicsTrigger() )
+ fNev++;
+ break;
+ case kCountAll:
+ fNev++;
+ break;
+ case kCountRaw:
+ fNev = evnum;
+ break;
+ default:
+ break;
+ }
+
+ //--- Print marks periodically
+ if( (fVerbose>1) && (evnum > 0) && (evnum % fMarkInterval == 0)){
+ cout << "test run: "<< fRun->GetNumber() << ", " << dec << evnum << endl;
+ }
+ if( (evnum > 0) &&(evnum % fAutoSaveInterval == 0)){
+ fOutput->GetTree()->AutoSave("SaveSelf");
+ }
+
+ //--- Update run parameters with current event
+ if( fUpdateRun ){
+ fRun->Update( fEvData );
+ }
+
+ //--- Clear all tests/cuts
+ if( fDoBench ) fBench->Begin("Cuts");
+ gHaCuts->ClearAll();
+ if( fDoBench ) fBench->Stop("Cuts");
+
+ //--- Perform the analysis
+ Int_t err = MainAnalysis();
+ switch( err ) {
+ case kOK:
+ break;
+ case kSkip:
+ continue;
+ case kFatal:
+ fatal = terminate = true;
+ continue;
+ case kTerminate:
+ terminate = true;
+ break;
+ default:
+ Error( here, "Unknown return code from MainAnalysis(): %d", err );
+ terminate = fatal = true;
+ continue;
+ }
+
+ Incr(kNevAccepted);
+
+ } // End of event loop
+
+ signal (SIGINT, prev_handler);
+ EndAnalysis();
+
+ //--- Close the input file
+ fRun->Close();
+
+ // Save final run parameters in run object of caller, if any
+ *run = *fRun;
+
+ // Write the output file and clean up.
+ // This writes the Tree as well as any objects (histograms etc.)
+ // that are defined in the current directory.
+
+ if( fDoBench ) fBench->Begin("Output");
+ // Ensure that we are in the output file's current directory
+ // ... someone might have pulled the rug from under our feet
+
+ // get the CURRENT file, since splitting might have occurred
+ if( fOutput && fOutput->GetTree() )
+ fFile = fOutput->GetTree()->GetCurrentFile();
+ if( fFile ) fFile->cd();
+ if( fOutput ) fOutput->End();
+ if( fFile ) {
+ fRun->Write("Run_Data"); // Save run data to ROOT file
+ // fFile->Write();//already done by fOutput->End()
+ fFile->Purge(); // get rid of excess object "cycles"
+ }
+ if( fDoBench ) fBench->Stop("Output");
+
+ fBench->Stop("Total");
+
+ //--- Report statistics
+ if( fVerbose>0 ) {
+ cout << dec;
+ if( status == THaRunBase::READ_EOF )
+ cout << "End of file";
+ else if ( fNev == nlast )
+ cout << "Event limit reached.";
+ else if ( fatal )
+ cout << "Fatal processing error.";
+ else if ( terminate )
+ cout << "Terminated during processing.";
+ cout << endl;
+
+ if( !fatal ) {
+ PrintCounters();
+
+ if( fVerbose>1 )
+ PrintScalers();
+ }
+ }
+
+ // Print cut summary (also to file if one given)
+ //if( !fatal )
+ // PrintCutSummary();
+
+ // Print timing statistics, if benchmarking enabled
+ //if( fDoBench && !fatal ) {
+ // cout << "Timing summary:" << endl;
+ // fBench->Print("Init");
+ // fBench->Print("RawDecode");
+ // fBench->Print("Decode");
+ // fBench->Print("CoarseTracking");
+ // fBench->Print("CoarseReconstruct");
+ // fBench->Print("Tracking");
+ // fBench->Print("Reconstruct");
+ // fBench->Print("Physics");
+ // fBench->Print("Output");
+ // fBench->Print("Cuts");
+ //}
+ //if( (fVerbose>1 || fDoBench) && !fatal )
+ // fBench->Print("Total");
+
+ //keep the last run available
+ // gHaRun = NULL;
+ return fNev;
+}
+
+}
diff --git a/src/Scandalizer.h b/src/Scandalizer.h
new file mode 100644
index 0000000000000000000000000000000000000000..7b08e1c1390f0220b46ab670f8a159a3a36759e5
--- /dev/null
+++ b/src/Scandalizer.h
@@ -0,0 +1,24 @@
+#ifndef hcana_Scandalizer_h_
+#define hcana_Scandalizer_h_
+
+#include "THaBenchmark.h"
+#include "THcAnalyzer.h"
+#include <iostream>
+
+namespace hcana {
+
+ class Scandalizer : public THcAnalyzer {
+ public:
+ Scandalizer() : THcAnalyzer() {}
+ virtual ~Scandalizer() {}
+
+ virtual Int_t Process(THaRunBase* run = nullptr);
+ virtual Int_t ReadOneEvent();
+
+ Int_t GoToEndOfCodaFile();
+
+ ClassDef(Scandalizer, 0) // Hall C Analyzer Standard Event Loop
+ };
+
+} // namespace hcana
+#endif
diff --git a/src/THcAnalyzer.h b/src/THcAnalyzer.h
index d5bd10bae17f11b4980e128a8d0745ae4c8f335b..80208f02d1b44ed4153f98f7c0dee446c80447e6 100644
--- a/src/THcAnalyzer.h
+++ b/src/THcAnalyzer.h
@@ -8,6 +8,41 @@
//////////////////////////////////////////////////////////////////////////
#include "THaAnalyzer.h"
+#include "THaRunBase.h"
+#include "THaEvent.h"
+#include "THaOutput.h"
+#include "THaEvData.h"
+#include "THaGlobals.h"
+#include "THaSpectrometer.h"
+#include "THaNamedList.h"
+#include "THaCutList.h"
+#include "THaCut.h"
+#include "THaPhysicsModule.h"
+#include "THaPostProcess.h"
+#include "THaBenchmark.h"
+#include "THaEvtTypeHandler.h"
+#include "THaEpicsEvtHandler.h"
+#include "TList.h"
+#include "TTree.h"
+#include "TFile.h"
+#include "TClass.h"
+#include "TDatime.h"
+#include "TClass.h"
+#include "TError.h"
+#include "TSystem.h"
+#include "TROOT.h"
+#include "TMath.h"
+#include "TDirectory.h"
+#include "THaCrateMap.h"
+
+#include <algorithm>
+#include <csignal>
+#include <cstring>
+#include <exception>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <stdexcept>
class THcAnalyzer : public THaAnalyzer {
diff --git a/src/THcHitList.cxx b/src/THcHitList.cxx
index 8422139db38f86c6dc2f2e238db95cc394349346..dce91765b5d04c85f9d41ecaf4ce85ae2c27ea16 100644
--- a/src/THcHitList.cxx
+++ b/src/THcHitList.cxx
@@ -186,7 +186,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
if(!fMap) { // Find the TI slot for ADCs
// Assumes that all FADCs are in the same crate
- cout << "Got the Crate map" << endl;
+ //cout << "Got the Crate map" << endl;
fMap = evdata.GetCrateMap();
for (Int_t i=0; i < fdMap->GetSize(); i++) { // Look for a FADC250
THaDetMap::Module* d = fdMap->GetModule(i);
@@ -197,7 +197,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
if(fMap->getModel(d->crate, slot) == 4) {
fTISlot = slot;
fTICrate = d->crate;
- cout << "TI Slot = " << fTISlot << endl;
+ //cout << "TI Slot = " << fTISlot << endl;
break;
}
}
diff --git a/src/TrackingEfficiency.cxx b/src/TrackingEfficiency.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..53b055a53a4d19b9c5326e8f34ec3aff765292b0
--- /dev/null
+++ b/src/TrackingEfficiency.cxx
@@ -0,0 +1,451 @@
+#include "THaEvData.h"
+#include "THaCutList.h"
+#include "VarDef.h"
+#include "VarType.h"
+#include "TClonesArray.h"
+
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <iostream>
+
+#include "THaApparatus.h"
+#include "THcGlobals.h"
+#include "THcHodoHit.h"
+#include "THcParmList.h"
+#include "TrackingEfficiency.h"
+
+namespace hcana {
+
+ using namespace std;
+
+ TrackingEfficiency::TrackingEfficiency(const char* name, const char* description,
+ const char* hodname)
+ : THaPhysicsModule(name, description) {}
+
+ TrackingEfficiency::~TrackingEfficiency() {
+ // Destructor
+ RemoveVariables();
+ }
+ //_____________________________________________________________________________
+
+ void TrackingEfficiency::Reset(Option_t* opt)
+ // Clear event-by-event data
+ {
+ Clear(opt);
+ }
+
+ //_____________________________________________________________________________
+ Int_t TrackingEfficiency::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;
+ // fBothGood[ip][ic] = 0;
+ // fPosGood[ip][ic] = 0;
+ // fNegGood[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 TrackingEfficiency::End(THaRunBase*) {
+ //// End of analysis
+ //for (Int_t ip = 0; ip < fNPlanes; ip++) {
+ // fStatAndEff[ip] = 0;
+ // for (Int_t ic = 0; ic < fNCounters[ip]; ic++) {
+ // fStatTrkSum[ip] += fStatTrk[fHod->GetScinIndex(ip, ic)];
+ // fStatAndSum[ip] += fHodoAndEffi[fHod->GetScinIndex(ip, ic)];
+ // }
+ // if (fStatTrkSum[ip] != 0)
+ // fStatAndEff[ip] = float(fStatAndSum[ip]) / float(fStatTrkSum[ip]);
+ //}
+ ////
+ //Double_t p1 = fStatAndEff[0];
+ //Double_t p2 = fStatAndEff[1];
+ //Double_t p3 = fStatAndEff[2];
+ //Double_t p4 = fStatAndEff[3];
+ //// probability that ONLY the listed planes had triggers
+ //Double_t p1234 = p1 * p2 * p3 * p4;
+ //Double_t p123 = p1 * p2 * p3 * (1. - p4);
+ //Double_t p124 = p1 * p2 * (1. - p3) * p4;
+ //Double_t p134 = p1 * (1. - p2) * p3 * p4;
+ //Double_t p234 = (1. - p1) * p2 * p3 * p4;
+ //fHodoEff_s1 = 1. - ((1. - p1) * (1. - p2));
+ //fHodoEff_s2 = 1. - ((1. - p3) * (1. - p4));
+ //fHodoEff_tof = fHodoEff_s1 * fHodoEff_s2;
+ //fHodoEff_3_of_4 = p1234 + p123 + p124 + p134 + p234;
+ //fHodoEff_4_of_4 = p1234;
+ return 0;
+ }
+
+ //_____________________________________________________________________________
+ THaAnalysisObject::EStatus TrackingEfficiency::Init(const TDatime& run_time) {
+ // Initialize TrackingEfficiency 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 << "TrackingEfficiency::Init nplanes=" << fHod->GetNPlanes() << endl;
+ //cout << "TrackingEfficiency::Init Apparatus = " << fHod->GetName() << " "
+ // << (fHod->GetApparatus())->GetName() << endl;
+
+ return fStatus = kOK;
+ }
+
+ //_____________________________________________________________________________
+ Int_t TrackingEfficiency::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];
+ //fStatTrkSum = new Int_t[fNPlanes];
+ //fStatAndSum = new Int_t[fNPlanes];
+ //fStatAndEff = new Double_t[fNPlanes];
+
+ //Int_t maxcountersperplane = 0;
+ //for (Int_t ip = 0; ip < fNPlanes; ip++) {
+ // fStatTrkSum[ip] = 0.;
+ // fStatAndSum[ip] = 0.;
+ // fStatAndEff[ip] = 0.;
+ // 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) + fPlanes[ip]->GetPosOffset();
+ // fNCounters[ip] = fPlanes[ip]->GetNelem();
+ // maxcountersperplane = TMath::Max(maxcountersperplane, fNCounters[ip]);
+ //}
+ //Int_t totalpaddles = fNPlanes * maxcountersperplane;
+ //fHodoPosEffi = new Int_t[totalpaddles];
+ //fHodoNegEffi = new Int_t[totalpaddles];
+ //fHodoOrEffi = new Int_t[totalpaddles];
+ //fHodoAndEffi = new Int_t[totalpaddles];
+ //fStatTrk = new Int_t[totalpaddles];
+
+ //char prefix[2];
+ //prefix[0] = tolower((fHod->GetApparatus())->GetName()[0]);
+ //prefix[1] = '\0';
+
+ //DBRequest list[] = {{"stat_slop", &fStatSlop, kDouble},
+ // {"stat_maxchisq", &fMaxChisq, kDouble},
+ // {"HodoEff_CalEnergy_Cut", &fHodoEff_CalEnergy_Cut, kDouble, 0, 1},
+ // {"hodo_slop", fHodoSlop, kDouble, (UInt_t)fNPlanes},
+ // {0}};
+ //fHodoEff_CalEnergy_Cut = 0.050; // set default value
+ //gHcParms->LoadParmValues((DBRequest*)&list, prefix);
+ //cout << "\n\nTrackingEfficiency::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);
+ //fBothGood.resize(fNPlanes);
+ //fPosGood.resize(fNPlanes);
+ //fNegGood.resize(fNPlanes);
+
+ //for (Int_t ip = 0; ip < fNPlanes; ip++) {
+
+ // cout << "Plane = " << ip + 1 << " counters = " << fNCounters[ip] << endl;
+
+ // 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]);
+ // fBothGood[ip].resize(fNCounters[ip]);
+ // fPosGood[ip].resize(fNCounters[ip]);
+ // fNegGood[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
+
+ // fHodoPosEffi[fHod->GetScinIndex(ip, ic)] = 0;
+ // fHodoNegEffi[fHod->GetScinIndex(ip, ic)] = 0;
+ // fHodoOrEffi[fHod->GetScinIndex(ip, ic)] = 0;
+ // fHodoAndEffi[fHod->GetScinIndex(ip, ic)] = 0;
+ // fStatTrk[fHod->GetScinIndex(ip, ic)] = 0;
+ // }
+ //}
+
+ //// Int_t fHodPaddles = fNCounters[0];
+ //// gHcParms->Define(Form("%shodo_pos_hits[%d][%d]",fPrefix,fNPlanes,fHodPaddles),
+ //// "Golden track's pos pmt hit",*&fStatPosHit);
+
+ //gHcParms->Define(Form("%shodo_pos_eff[%d]", prefix, totalpaddles), "Hodo positive effi",
+ // *fHodoPosEffi);
+ //gHcParms->Define(Form("%shodo_neg_eff[%d]", prefix, totalpaddles), "Hodo negative effi",
+ // *fHodoNegEffi);
+ //gHcParms->Define(Form("%shodo_or_eff[%d]", prefix, totalpaddles), "Hodo or effi", *fHodoOrEffi);
+ //gHcParms->Define(Form("%shodo_and_eff[%d]", prefix, totalpaddles), "Hodo and effi",
+ // *fHodoAndEffi);
+ //gHcParms->Define(Form("%shodo_plane_AND_eff[%d]", prefix, fNPlanes), "Hodo plane AND eff",
+ // *fStatAndEff);
+ //gHcParms->Define(Form("%shodo_gold_hits[%d]", prefix, totalpaddles), "Hodo golden hits",
+ // *fStatTrk);
+ //gHcParms->Define(Form("%shodo_s1XY_eff", prefix), "Efficiency for S1XY", fHodoEff_s1);
+ //gHcParms->Define(Form("%shodo_s2XY_eff", prefix), "Efficiency for S2XY", fHodoEff_s2);
+ //gHcParms->Define(Form("%shodo_stof_eff", prefix), "Efficiency for STOF", fHodoEff_tof);
+ //gHcParms->Define(Form("%shodo_3_of_4_eff", prefix), "Efficiency for 3 of 4", fHodoEff_3_of_4);
+ //gHcParms->Define(Form("%shodo_4_of_4_eff", prefix), "Efficiency for 4 of 4", fHodoEff_4_of_4);
+
+ return kOK;
+ }
+
+ //_____________________________________________________________________________
+ Int_t TrackingEfficiency::DefineVariables(EMode mode) {
+
+ //if (mode == kDefine && fIsSetup)
+ // return kOK;
+ //fIsSetup = (mode == kDefine);
+
+ //// fEffiTest = 0;
+ //// gHcParms->Define(Form("hodoeffi"),"Testing effi",fEffiTest);
+
+ //const RVarDef vars[] = {// Move these into THcHallCSpectrometer using track fTracks
+ // // {"effitestvar", "efficiency test var", "fEffiTest"},
+ // // {"goldhodposhit", "pos pmt hit in hodo", "fStatPosHit"},
+ // {0}};
+ //return DefineVarsFromList(vars, mode);
+ return kOK;
+ }
+
+ //_____________________________________________________________________________
+ Int_t TrackingEfficiency::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();
+ //// 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;
+ //Int_t trackIndex = theTrack->GetTrkNum() - 1;
+
+ //// 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];
+ // hitCounter[ip] =
+ // TMath::Max(TMath::Min(TMath::Nint((hitPos[ip] - fCenterFirst[ip]) / fSpacing[ip] + 1),
+ // fNCounters[ip]),
+ // 1);
+ // hitDistance[ip] = hitPos[ip] - (fSpacing[ip] * (hitCounter[ip] - 1) + fCenterFirst[ip]);
+ // } else { // Y Plane
+ // hitPos[ip] = theTrack->GetY() + theTrack->GetPhi() * fPosZ[ip];
+ // hitCounter[ip] =
+ // TMath::Max(TMath::Min(TMath::Nint((fCenterFirst[ip] - hitPos[ip]) / fSpacing[ip] + 1),
+ // fNCounters[ip]),
+ // 1);
+ // hitDistance[ip] = hitPos[ip] - (fCenterFirst[ip] - fSpacing[ip] * (hitCounter[ip] - 1));
+ // }
+ //}
+
+ //// 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 &&
+ // theTrack->GetEnergy() >= fHodoEff_CalEnergy_Cut) {
+ // fStatTrk[fHod->GetScinIndex(ip, hitCounter[ip] - 1)]++;
+ // // 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];
+ // if (hitcounter >= fNCounters[ip])
+ // hitcounter = fNCounters[ip] - 1;
+ // if (hitcounter < 0)
+ // hitcounter = 0;
+ // 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 &&
+ // theTrack->GetEnergy() >= fHodoEff_CalEnergy_Cut) {
+ // 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]++;
+
+ // fHodoPosEffi[fHod->GetScinIndex(ip, hitCounter[ip] - 1)]++;
+ // fHodoNegEffi[fHod->GetScinIndex(ip, hitCounter[ip] - 1)]++;
+ // fHodoAndEffi[fHod->GetScinIndex(ip, hitCounter[ip] - 1)]++;
+ // fHodoOrEffi[fHod->GetScinIndex(ip, hitCounter[ip] - 1)]++;
+
+ // // 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]++;
+ // fHodoPosEffi[fHod->GetScinIndex(ip, hitCounter[ip] - 1)]++;
+ // fHodoOrEffi[fHod->GetScinIndex(ip, hitCounter[ip] - 1)]++;
+ // }
+ // } else if (goodTdcNeg) {
+ // fStatNegHit[ip][hitcounter]++;
+ // fStatOrHit[ip][hitcounter]++;
+ // fHodoNegEffi[fHod->GetScinIndex(ip, hitCounter[ip] - 1)]++;
+ // fHodoOrEffi[fHod->GetScinIndex(ip, hitCounter[ip] - 1)]++;
+ // }
+
+ // // 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;
+ }
+
+} // namespace hcana
+
diff --git a/src/TrackingEfficiency.h b/src/TrackingEfficiency.h
new file mode 100644
index 0000000000000000000000000000000000000000..5e2ea1b5866891e87ab41c6c8a119eaf9af941f3
--- /dev/null
+++ b/src/TrackingEfficiency.h
@@ -0,0 +1,104 @@
+#ifndef ROOT_TrackingEfficiency
+#define ROOT_TrackingEfficiency
+
+///////////////////////////////////////////////////////////////////////////////
+// //
+// TrackingEfficiency //
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "THaEvData.h"
+#include "THaCutList.h"
+#include "VarDef.h"
+#include "VarType.h"
+#include "TClonesArray.h"
+
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+
+#include "THaPhysicsModule.h"
+#include "THcHodoscope.h"
+#include "THaSpectrometer.h"
+#include "THaTrack.h"
+
+namespace hcana {
+
+
+ /** \brief TrackingEfficiency calculation.
+ *
+ * \ingroup PhysMods
+ */
+ class TrackingEfficiency : public THaPhysicsModule {
+ public:
+ TrackingEfficiency(const char* name, const char* description, const char* hodname);
+ virtual ~TrackingEfficiency();
+
+ 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);
+ virtual Int_t DefineVariables(EMode mode = kDefine);
+ /* Int_t GetScinIndex(Int_t nPlane, Int_t nPaddle); */
+
+ // 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 fEffiTest;
+ // Int_t fNPlanes;
+ // THcScintillatorPlane** fPlanes;
+ // Double_t* fPosZ;
+ // Double_t* fSpacing;
+ // Double_t* fCenterFirst;
+ // Int_t* fNCounters;
+ //// Int_t* fHodoPlnContHit;
+ // Int_t* fHodoPosEffi;
+ // Int_t* fHodoNegEffi;
+ // Int_t* fHodoOrEffi;
+ // Int_t* fHodoAndEffi;
+ // Int_t* fStatTrk;
+ // Int_t* fStatTrkSum;
+ // Int_t* fStatAndSum;
+ // Double_t* fStatAndEff;
+ // Double_t fStatSlop;
+ // Double_t fHodoEff_CalEnergy_Cut;
+ // Double_t fMaxChisq;
+ // Double_t* fHodoSlop;
+ // Double_t fHodoEff_s1, fHodoEff_s2, fHodoEff_tof, fHodoEff_3_of_4, fHodoEff_4_of_4;
+
+ //// 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(TrackingEfficiency, 0) // Hodoscope efficiency module
+ };
+
+} // namespace hcana
+
+#endif
diff --git a/src/include/HallC_LinkDef.h b/src/include/HallC_LinkDef.h
index d3af17de665afe8abbb0f3355f7bf96a89c98c3f..e666b7e61449aa49fb75adc7f697ceeecd9b05f4 100644
--- a/src/include/HallC_LinkDef.h
+++ b/src/include/HallC_LinkDef.h
@@ -77,6 +77,8 @@
#pragma link C++ class Decoder::TIBlobModule+;
#pragma link C++ class hcana::Shower2+;
+#pragma link C++ class hcana::Scandalizer+;
+#pragma link C++ class hcana::TrackingEfficiency+;
// Postamble for HallC_Linkdef.h file
#endif