THcHodoscope.cxx

/** \class THcHodoscope
    \ingroup Detectors

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 "THcSignalHit.h"
#include "THcShower.h"
#include "THcCherenkov.h"
#include "THcHallCSpectrometer.h"

#include "THcHitList.h"
#include "THcRawShowerHit.h"
#include "TClass.h"
#include "math.h"
#include "THaSubDetector.h"

#include "THcHodoscope.h"
#include "THaEvData.h"
#include "THaDetMap.h"
#include "THcDetectorMap.h"
#include "THaGlobals.h"
#include "THaCutList.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 <vector>

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

using namespace std;
using std::vector;

//_____________________________________________________________________________
THcHodoscope::THcHodoscope( const char* name, const char* description,
				  THaApparatus* apparatus ) :
  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( ) :
  THaNonTrackingDetector()
{
  // Constructor
}

//_____________________________________________________________________________
void THcHodoscope::Setup(const char* name, const char* description)
{

  //  static const char* const here = "Setup()";
  //  static const char* const message =
  //    "Must construct %s detector with valid name! Object construction failed.";

  cout << "In THcHodoscope::Setup()" << endl;
  // Base class constructor failed?
  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;
  cout << " fSHMS = " << fSHMS << endl;
  string planenamelist;
  DBRequest listextra[]={
    {"hodo_num_planes", &fNPlanes, kInt},
    {"hodo_plane_names",&planenamelist, kString},
    {0}
  };
  //fNPlanes = 4; 		// Default if not defined
  gHcParms->LoadParmValues((DBRequest*)&listextra,prefix);

  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;
  }

  // --------------- To get energy from THcShower ----------------------
  const char* shower_detector_name = "cal";
  //  THaApparatus* app;
  THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
  THaDetector* det = app->GetDetector( shower_detector_name );

  if( dynamic_cast<THcShower*>(det) ) {
    fShower = dynamic_cast<THcShower*>(det);
  }
  else if( !dynamic_cast<THcShower*>(det) ) {
    cout << "Warining: calorimeter analysis module "
	 << shower_detector_name << " not loaded for spectrometer "
	 << prefix << endl;

    fShower = NULL;
  }

  // --------------- To get energy from THcShower ----------------------

  // --------------- To get NPEs from THcCherenkov -------------------
  const char* chern_detector_name = "cer";
  THaDetector* detc = app->GetDetector( chern_detector_name );

  if( dynamic_cast<THcCherenkov*>(detc) ) {
    fChern = dynamic_cast<THcCherenkov*>(detc);
  }
  else if( !dynamic_cast<THcCherenkov*>(detc) ) {
    cout << "Warining: Cherenkov detector analysis module "
	 << chern_detector_name << " not loaded for spectrometer "
	 << prefix << endl;

    fChern = NULL;
  }

  // --------------- To get NPEs from THcCherenkov -------------------

  fScinShould = 0;
  fScinDid = 0;
  gHcParms->Define(Form("%shodo_did",prefix),"Total hodo tracks",fScinDid);
  gHcParms->Define(Form("%shodo_should",prefix),"Total hodo triggers",fScinShould);

  // Save the nominal particle mass
  fPartMass = app->GetParticleMass();
  fBetaNominal = app->GetBetaAtPcentral();

  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 );
    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.

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

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