THcShower.cxx

/** \class THcShower
    \ingroup Detectors

Shower counter class, describing a generic segmented shower detector.     //

*/

#include "THcShower.h"
#include "THcHallCSpectrometer.h"
#include "THaEvData.h"
#include "THaDetMap.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
#include "THaCutList.h"
#include "THcParmList.h"
#include "VarDef.h"
#include "VarType.h"
#include "THaTrack.h"
#include "TClonesArray.h"
#include "THaTrackProj.h"
#include "TMath.h"

#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <numeric>

using namespace std;

//_____________________________________________________________________________
THcShower::THcShower( const char* name, const char* description,
				  THaApparatus* apparatus ) :
  THaNonTrackingDetector(name,description,apparatus)
{
  // Constructor
  fNLayers = 0;			// No layers until we make them
  fNTotLayers = 0;
  fHasArray = 0;

  fClusterList = new THcShowerClusterList;
}

//_____________________________________________________________________________
THcShower::THcShower( ) :
  THaNonTrackingDetector()
{
  // Constructor
}

//_____________________________________________________________________________
void THcShower::Setup(const char* name, const char* description)
{
  char prefix[2];

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

  
  string layernamelist;
  fHasArray = 0;		// Flag for presence of fly's eye array
  DBRequest list[]={
    {"cal_num_layers", &fNLayers, kInt},
    {"cal_layer_names", &layernamelist, kString},
    {"cal_array",&fHasArray, kInt,0, 1}, 
    {0}
  };

  gHcParms->LoadParmValues((DBRequest*)&list,prefix);
  fNTotLayers = (fNLayers+(fHasArray!=0?1:0));
  vector<string> layer_names = vsplit(layernamelist);

  if(layer_names.size() != fNTotLayers) {
    cout << "THcShower::Setup ERROR: Number of layers " << fNTotLayers 
	 << " doesn't agree with number of layer names "
	 << layer_names.size() << endl;
    // Should quit.  Is there an official way to quit?
  }

  fLayerNames = new char* [fNTotLayers];
  for(UInt_t i=0;i<fNTotLayers;i++) {
    fLayerNames[i] = new char[layer_names[i].length()+1];
    strcpy(fLayerNames[i], layer_names[i].c_str());
  }
  
  char *desc = new char[strlen(description)+100];
  fPlanes = new THcShowerPlane* [fNLayers];

  for(UInt_t i=0;i < fNLayers;i++) {
    strcpy(desc, description);
    strcat(desc, " Plane ");
    strcat(desc, fLayerNames[i]);

    fPlanes[i] = new THcShowerPlane(fLayerNames[i], desc, i+1, this); 

  }
  if(fHasArray) {
    strcpy(desc, description);
    strcat(desc, " Array ");
    strcat(desc, fLayerNames[fNTotLayers-1]);

    fArray = new THcShowerArray(fLayerNames[fNTotLayers-1], desc, fNTotLayers,
				this);
  } else {
    fArray = 0;
  }

  delete [] desc;

  cout << "---------------------------------------------------------------\n";

  cout << "From THcShower::Setup: created Shower planes for "
       << GetApparatus()->GetName() << ": ";

  for(UInt_t i=0;i < fNTotLayers;i++) {
    cout << fLayerNames[i];
    i < fNTotLayers-1 ? cout << ", " : cout << ".\n";
  }

  if(fHasArray)
    cout << fLayerNames[fNTotLayers-1] << " has fly\'s eye configuration\n";

  cout << "---------------------------------------------------------------\n";

}


//_____________________________________________________________________________
THaAnalysisObject::EStatus THcShower::Init( const TDatime& date )
{
  Setup(GetName(), GetTitle());

  // Should probably put this in ReadDatabase as we will know the
  // maximum number of hits after setting up the detector map

  InitHitList(fDetMap, "THcRawShowerHit", 100);

  EStatus status;
  if( (status = THaNonTrackingDetector::Init( date )) )
    return fStatus=status;

  for(UInt_t ip=0;ip<fNLayers;ip++) {
    if((status = fPlanes[ip]->Init( date ))) {
      return fStatus=status;
    }
  }
  if(fHasArray) {
    if((status = fArray->Init( date ))) {
      return fStatus = status;
    }
  }

  char EngineDID[] = " CAL";
  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;
  }

  if(fHasArray) {
    cout << "THcShower::Init: adjustment of fiducial volume limits to the fly's eye part." << endl;
    cout << "  Old limits:" << endl;
    cout << "    Xmin = " << fvXmin << "  Xmax = " << fvXmax << endl;
    cout << "    Ymin = " << fvYmin << "  Ymax = " << fvYmax << endl;
    fvXmin = TMath::Max(fvXmin, fArray->fvXmin());
    fvXmax = TMath::Min(fvXmax, fArray->fvXmax());
    fvYmin = TMath::Max(fvYmin, fArray->fvYmin());
    fvYmax = TMath::Min(fvYmax, fArray->fvYmax());
    cout << "  New limits:" << endl;
    cout << "    Xmin = " << fvXmin << "  Xmax = " << fvXmax << endl;
    cout << "    Ymin = " << fvYmin << "  Ymax = " << fvYmax << endl;
  }

  cout << "---------------------------------------------------------------\n";
  cout << "From THcShower::Init: initialized " << GetApparatus()->GetName()
       <<  GetName() << endl;
  cout << "---------------------------------------------------------------\n";

  return fStatus = kOK;
}

//_____________________________________________________________________________
Int_t THcShower::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];

  // 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