THcHodoscope.cxx

/** \class THcHodoscope
    \ingroup Detectors

\brief 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 "TClass.h"
#include "THaSubDetector.h"
#include "THcCherenkov.h"
#include "THcHallCSpectrometer.h"
#include "THcHitList.h"
#include "THcRawShowerHit.h"
#include "THcScintPlaneCluster.h"
#include "THcShower.h"
#include "THcSignalHit.h"
#include "math.h"

#include "TClonesArray.h"
#include "THaCutList.h"
#include "THaDetMap.h"
#include "THaEvData.h"
#include "THaGlobals.h"
#include "THaTrack.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
#include "THcHodoscope.h"
#include "THcParmList.h"
#include "TMath.h"
#include "VarDef.h"
#include "VarType.h"

#include "THaOutput.h"
#include "TTree.h"

#include "THaTrackProj.h"
#include <vector>

#include <algorithm>
#include <array>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>

#include "hcana/helpers.hxx"
#include <cassert>

using namespace std;
using std::vector;

//_____________________________________________________________________________
THcHodoscope::THcHodoscope(const char* name, const char* description, THaApparatus* apparatus)
    : hcana::ConfigLogging<THaNonTrackingDetector>(name, description, apparatus) {
  // Constructor

  // fTrackProj = new TClonesArray( "THaTrackProj", 5 );
  // Construct the planes
  fNPlanes       = 0; // No planes until we make them
  fStartTime     = -1e5;
  fGoodStartTime = kFALSE;
}

//_____________________________________________________________________________
THcHodoscope::THcHodoscope() : hcana::ConfigLogging<THaNonTrackingDetector>() {
  // Constructor
}

//_____________________________________________________________________________
void THcHodoscope::Setup(const char* name, const char* description) {
  /**
     Create the scintillator plane objects for the hodoscope.

     Uses the Xhodo_num_planes and Xhodo_plane_names to get the number of
     planes and their names.

     Gets a pointer to the Cherenkov named "cer" ("hgcer" in the case of the SHMS.)

  */
  if (IsZombie())
    return;

  // fDebug = 1;  // Keep this at one while we're working on the code

  char prefix[2];

  prefix[0] = tolower(GetApparatus()->GetName()[0]);
  prefix[1] = '\0';

  TString temp(prefix[0]);
  fSHMS = kFALSE;
  if (temp == "p")
    fSHMS = kTRUE;
  TString histname = temp + "_timehist";
  hTime            = new TH1F(histname, "", 400, 0, 200);
  // cout << " fSHMS = " << fSHMS << endl;
  string    planenamelist;
  DBRequest listextra[] = {{"hodo_num_planes", &fNPlanes, kInt},
                           {"hodo_plane_names", &planenamelist, kString},
                           {"hodo_tdcrefcut", &fTDC_RefTimeCut, kInt, 0, 1},
                           {"hodo_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
                           {0}};
  // fNPlanes = 4; 		// Default if not defined
  fTDC_RefTimeCut = 0; // Minimum allowed reference times
  fADC_RefTimeCut = 0;
  gHcParms->LoadParmValues((DBRequest*)&listextra, prefix);

  _det_logger->info("Plane Name List : {}", planenamelist);
  // cout << "Plane Name List : " << planenamelist << endl;

  vector<string> plane_names = vsplit(planenamelist);
  // Plane names
  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() + 1];
    strcpy(fPlaneNames[i], plane_names[i].c_str());
  }

  // Probably shouldn't assume that description is defined
  char* desc = new char[strlen(description) + 100];
  fPlanes    = new THcScintillatorPlane*[fNPlanes];
  for (Int_t i = 0; i < fNPlanes; i++) {
    strcpy(desc, description);
    strcat(desc, " Plane ");
    strcat(desc, fPlaneNames[i]);
    fPlanes[i] = new THcScintillatorPlane(fPlaneNames[i], desc, i + 1,
                                          this); // Number planes starting from zero!!
    // cout << "Created Scintillator Plane " << fPlaneNames[i] << ", " << desc << endl;
  }

  // Save the nominal particle mass
  THcHallCSpectrometer* app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
  fPartMass                 = app->GetParticleMass();
  fBetaNominal              = app->GetBetaAtPcentral();

  if (fSHMS) {
    fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("hgcer"));
  } else {
    fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("cer"));
  }

  delete[] desc;
}

//_____________________________________________________________________________
THaAnalysisObject::EStatus THcHodoscope::Init(const TDatime& date) {
  // cout << "In THcHodoscope::Init()" << endl;
  Setup(GetName(), GetTitle());

  char EngineDID[] = "xSCIN";
  EngineDID[0]     = toupper(GetApparatus()->GetName()[0]);
  if (gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0) {
    static const char* const here = "Init()";
    // Error( Here(here), "Error filling detectormap for %s.", EngineDID );
    _det_logger->error("THcHodoscope::Init : Error filling detectormap for {}.", EngineDID);
    return kInitError;
  }

  // Should probably put this in ReadDatabase as we will know the
  // maximum number of hits after setting up the detector map
  // But it needs to happen before the sub detectors are initialized
  // so that they can get the pointer to the hitlist.
  _det_logger->info("Hodo TDC and ADC ref time cut = {} {}", fTDC_RefTimeCut, fADC_RefTimeCut);
  // cout << " Hodo tdc ref time cut = " << fTDC_RefTimeCut << " " << fADC_RefTimeCut << endl;

  InitHitList(fDetMap, "THcRawHodoHit", fDetMap->GetTotNumChan() + 1, fTDC_RefTimeCut,
              fADC_RefTimeCut);

  EStatus status;
  // This triggers call of ReadDatabase and DefineVariables
  if ((status = THaNonTrackingDetector::Init(date)))
    return fStatus = status;

  for (Int_t ip = 0; ip < fNPlanes; ip++) {
    if ((status = fPlanes[ip]->Init(date))) {
      return fStatus = status;
    }
  }

  // Why not std::vector?
  fNScinHits     = new Int_t[fNPlanes];
  fGoodPlaneTime = new Bool_t[fNPlanes];
  fNPlaneTime    = new Int_t[fNPlanes];
  fSumPlaneTime  = new Double_t[fNPlanes];

  //  Double_t  fHitCnt4 = 0., fHitCnt3 = 0.;

  // Int_t m = 0;
  // fScinHit = new Double_t*[fNPlanes];