THcShower.cxx

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

#include "THcShower.h"
#include "THaEvData.h"
#include "THaDetMap.h"
#include "THcDetectorMap.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 <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>

using namespace std;

//_____________________________________________________________________________
THcShower::THcShower( const char* name, const char* description,
				  THaApparatus* apparatus ) :
  THaNonTrackingDetector(name,description,apparatus)
{
  // Constructor

//  fTrackProj = new TClonesArray( "THaTrackProj", 5 );
}

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

//_____________________________________________________________________________
THaAnalysisObject::EStatus THcShower::Init( const TDatime& date )
{
  static const char* const here = "Init()";

  if( THaNonTrackingDetector::Init( date ) )
    return fStatus;


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

  THcHitList::InitHitList(fDetMap, "THcShowerHit", 100);

  // Will need to determine which apparatus it belongs to and use the
  // appropriate detector ID in the FillMap call
  if( gHcDetectorMap->FillMap(fDetMap, "HCAL") < 0 ) {
    Error( Here(here), "Error filling detectormap for %s.", 
	     "HCAL");
      return kInitError;
  }

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

  // Read data from database 
  // Pull values from the THcParmList instead of reading a database
  // file like Hall A does.

  //  DBRequest list[] = {
  //    { "TDC_offsetsL", fLOff, kDouble, fNelem },
  //    { "TDC_offsetsR", fROff, kDouble, fNelem },
  //    { "ADC_pedsL", fLPed, kDouble, fNelem },
  //    { "ADC_pedsR", fRPed, kDouble, fNelem },
  //    { "ADC_coefL", fLGain, kDouble, fNelem },
  //    { "ADC_coefR", fRGain, kDouble, fNelem },
  //    { "TDC_res",   &fTdc2T },
  //    { "TransSpd",  &fCn },
  //    { "AdcMIP",    &fAdcMIP },
  //    { "NTWalk",    &fNTWalkPar, kInt },
  //    { "Timewalk",  fTWalkPar, kDouble, 2*fNelem },
  //    { "ReTimeOff", fTrigOff, kDouble, fNelem },
  //    { "AvgRes",    &fResolution },
  //    { "Atten",     &fAttenuation },
  //    { 0 }
  //  };

  // We will probably want to add some kind of method to gHcParms to allow
  // bulk retrieval of parameters of interest.

  // Will need to determine which spectrometer in order to construct
  // the parameter names (e.g. hscin_1x_nr vs. sscin_1x_nr)

cout << "THcShower::ReadDatabase called " << GetName() << endl;

  NLayers = 4;			// Hardwire for now

  BlockThick = new Double_t [NLayers];

  BlockThick[0] = *(Double_t *)gHcParms->Find("hcal_1pr_thick")->GetValuePointer();
  BlockThick[1] = *(Double_t *)gHcParms->Find("hcal_2ta_thick")->GetValuePointer();
  BlockThick[2] = *(Double_t *)gHcParms->Find("hcal_3ta_thick")->GetValuePointer();
  BlockThick[3] = *(Double_t *)gHcParms->Find("hcal_4ta_thick")->GetValuePointer();

cout << "Block thickness: " << BlockThick[2]  << endl;

fNBlocks = new Int_t [NLayers];

  fNBlocks[0] = *(Int_t *)gHcParms->Find("hcal_1pr_nr")->GetValuePointer();
  fNBlocks[1] = *(Int_t *)gHcParms->Find("hcal_2ta_nr")->GetValuePointer();
  fNBlocks[2] = *(Int_t *)gHcParms->Find("hcal_3ta_nr")->GetValuePointer();
  fNBlocks[3] = *(Int_t *)gHcParms->Find("hcal_4ta_nr")->GetValuePointer();

  cout << "Number of blocks per layer: " << fNBlocks[2]  << endl;

  fNLayerZPos = new Double_t [NLayers];

  fNLayerZPos[0] = *(Double_t *)gHcParms->Find("hcal_1pr_zpos")->GetValuePointer();
  fNLayerZPos[1] = *(Double_t *)gHcParms->Find("hcal_2ta_zpos")->GetValuePointer();
  fNLayerZPos[2] = *(Double_t *)gHcParms->Find("hcal_3ta_zpos")->GetValuePointer();
  fNLayerZPos[3] = *(Double_t *)gHcParms->Find("hcal_4ta_zpos")->GetValuePointer();

cout << "Z Position: " << fNLayerZPos[2]  << endl;

  XPos = new Double_t [2*NLayers];

  XPos[0] = *(Double_t *)gHcParms->Find("hcal_1pr_left")->GetValuePointer();
  XPos[1] = *(Double_t *)gHcParms->Find("hcal_1pr_right")->GetValuePointer();
  XPos[2] = *(Double_t *)gHcParms->Find("hcal_2ta_left")->GetValuePointer();
  XPos[3] = *(Double_t *)gHcParms->Find("hcal_2ta_right")->GetValuePointer();
  XPos[4] = *(Double_t *)gHcParms->Find("hcal_3ta_left")->GetValuePointer();
  XPos[5] = *(Double_t *)gHcParms->Find("hcal_3ta_right")->GetValuePointer();
  XPos[6] = *(Double_t *)gHcParms->Find("hcal_4ta_left")->GetValuePointer();
  XPos[7] = *(Double_t *)gHcParms->Find("hcal_4ta_right")->GetValuePointer();

cout << "X Positions: " << XPos[0]  << ", " << XPos[1] << endl;

  YPos = new Double_t* [4];
  cout << "Y Positions:";

  Double_t* p;
  Int_t ilayer;

  ilayer = 0; 
  p = (Double_t *)gHcParms->Find("hcal_1pr_top")->GetValuePointer();
  YPos[ilayer] = new Double_t [fNBlocks[ilayer]];
  // Print out some parameters just to demonstrate that it works
  
  for(Int_t i=0;i<fNBlocks[ilayer];i++) {
    YPos[ilayer][i] = p[i];
    cout << " " << YPos[ilayer][i];
  }
  cout << endl;

  ilayer = 1; 
  p = (Double_t *)gHcParms->Find("hcal_2ta_top")->GetValuePointer();
  YPos[ilayer] = new Double_t [fNBlocks[ilayer]];
  
   for(Int_t i=0;i<fNBlocks[ilayer];i++) {
    YPos[ilayer][i] = p[i];
  }
 
  ilayer = 2; 
  p = (Double_t *)gHcParms->Find("hcal_3ta_top")->GetValuePointer();
  YPos[ilayer] = new Double_t [fNBlocks[ilayer]];

  for(Int_t i=0;i<fNBlocks[ilayer];i++) {
    YPos[ilayer][i] = p[i];
  }

  ilayer = 3; 
  p = (Double_t *)gHcParms->Find("hcal_4ta_top")->GetValuePointer();
  YPos[ilayer] = new Double_t [fNBlocks[ilayer]];

  for(Int_t i=0;i<fNBlocks[ilayer];i++) {
    YPos[ilayer][i] = p[i];
  }

  fIsInit = true;

  return kOK;
}



//_____________________________________________________________________________
Int_t THcShower::DefineVariables( EMode mode )
{
  // Initialize global variables and lookup table for decoder

  if( mode == kDefine && fIsSetup ) return kOK;
  fIsSetup = ( mode == kDefine );

  // Register variables in global list

 //  RVarDef vars[] = {
 //   { "nhit",   "Number of hits",                     "fNhits" },
 //   { "a",      "Raw ADC amplitudes",                 "fA" },
 //   { "a_p",    "Ped-subtracted ADC amplitudes",      "fA_p" },
 //   { "a_c",    "Calibrated ADC amplitudes",          "fA_c" },
 //   { "asum_p", "Sum of ped-subtracted ADCs",         "fAsum_p" },
 //   { "asum_c", "Sum of calibrated ADCs",             "fAsum_c" },
 //   { "nclust", "Number of clusters",                 "fNclust" },
 //   { "e",      "Energy (MeV) of largest cluster",    "fE" },
 //   { "x",      "x-position (cm) of largest cluster", "fX" },
 //   { "y",      "y-position (cm) of largest cluster", "fY" },
 //   { "mult",   "Multiplicity of largest cluster",    "fMult" },
 //   { "nblk",   "Numbers of blocks in main cluster",  "fNblk" },
 //   { "eblk",   "Energies of blocks in main cluster", "fEblk" },
 //   { "trx",    "track x-position in det plane",      "fTRX" },
 //   { "try",    "track y-position in det plane",      "fTRY" },
 //   { 0 }
 // };
 //  return DefineVarsFromList( vars, mode );
  return kOK;
}

//_____________________________________________________________________________
THcShower::~THcShower()
{
  // Destructor. Remove variables from global list.

  if( fIsSetup )
    RemoveVariables();
  if( fIsInit )
    DeleteArrays();
  if (fTrackProj) {
    fTrackProj->Clear();
    delete fTrackProj; fTrackProj = 0;
  }
}

//_____________________________________________________________________________
void THcShower::DeleteArrays()
{
  // Delete member arrays. Used by destructor.

  delete [] BlockThick;  BlockThick = NULL;
  delete [] fNBlocks;  fNBlocks = NULL;
  delete [] fNLayerZPos;  fNLayerZPos = NULL;
  delete [] XPos;  XPos = NULL;
  delete [] YPos;  YPos = NULL;
  //delete [] fSpacing;  fSpacing = NULL;
  //delete [] fCenter;   fCenter = NULL; // This 2D. What is correct way to delete?
}

//_____________________________________________________________________________
inline 
void THcShower::ClearEvent()
{
  // Reset per-event data.

  fTrackProj->Clear();
}

//_____________________________________________________________________________
Int_t THcShower::Decode( const THaEvData& evdata )
{

  // Get the Hall C style hitlist (fRawHitList) for this event
  Int_t nhits = THcHitList::DecodeToHitList(evdata);

//   fRawHitList is TClones array of THcShowerHit objects
  for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
    THcShowerHit* hit = (THcShowerHit *) fRawHitList->At(ihit);
    cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
	 << hit->fADC_pos << " " << hit->fADC_neg << " "  << endl;
  }
  cout << endl;

  return nhits;
}

//_____________________________________________________________________________
Int_t THcShower::ApplyCorrections( void )
{
  return(0);
}

//_____________________________________________________________________________
Double_t THcShower::TimeWalkCorrection(const Int_t& paddle,
					     const ESide side)
{
  return(0.0);
}

//_____________________________________________________________________________
Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
{
  // Calculation of coordinates of particle track cross point with shower
  // plane in the detector coordinate system. For this, parameters of track 
  // reconstructed in THaVDC::CoarseTrack() are used.
  //
  // Apply corrections and reconstruct the complete hits.
  //
  //  static const Double_t sqrt2 = TMath::Sqrt(2.);
  
  ApplyCorrections();

  return 0;
}

//_____________________________________________________________________________
Int_t THcShower::FineProcess( TClonesArray& tracks )
{
  // Reconstruct coordinates of particle track cross point with shower
  // plane, and copy the data into the following local data structure:
  //
  // Units of measurements are meters.

  // Calculation of coordinates of particle track cross point with shower
  // plane in the detector coordinate system. For this, parameters of track 
  // reconstructed in THaVDC::FineTrack() are used.

  return 0;
}

ClassImp(THcShower)
////////////////////////////////////////////////////////////////////////////////