THcShower.cxx

/** \class THcShower
    \ingroup Detectors

\brief 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},
    {"cal_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
    {0}
  };

  fADC_RefTimeCut = 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());

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

  InitHitList(fDetMap, "THcRawShowerHit", fDetMap->GetTotNumChan()+1,
	      0, fADC_RefTimeCut);

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

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

  fPresentP = 0;
  THaVar* vpresent = gHaVars->Find(Form("%s.present",GetApparatus()->GetName()));
  if(vpresent) {
    fPresentP = (Bool_t *) vpresent->GetValuePointer();
  }
  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];