-
Simon Zhamkochyan authoredSimon Zhamkochyan authored
THcShower.cxx 12.28 KiB
///////////////////////////////////////////////////////////////////////////////
// //
// 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
Setup(name, description);
// fTrackProj = new TClonesArray( "THaTrackProj", 5 );
}
//_____________________________________________________________________________
THcShower::THcShower( ) :
THaNonTrackingDetector()
{
// Constructor
}
void THcShower::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.";
// Base class constructor failed?
if( IsZombie()) return;
fNLayers = 4; // Eventually get # layers and layer names from a DB
fLayerNames = new char* [fNLayers];
for(Int_t i=0;i<fNLayers;i++) {fLayerNames[i] = new char[3];}
strcpy(fLayerNames[0],"1z");
strcpy(fLayerNames[1],"2z");
strcpy(fLayerNames[2],"3z");
strcpy(fLayerNames[3],"4z");
size_t nlen = strlen(name);
size_t slen = 0;
for(Int_t i=0;i < fNLayers;i++) {slen = TMath::Max(slen,strlen(fLayerNames[i]));}
size_t len = nlen+slen+1;
// Probably shouldn't assume that description is defined
char* desc = new char[strlen(description)+50+slen];
fPlanes = new THcShowerPlane* [fNLayers];
for(Int_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);
cout << "Created Shower Plane " << fLayerNames[i] << ", " << desc << endl;
}
}
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcShower::Init( const TDatime& date )
{
static const char* const here = "Init()";
cout << "THcShower::Init " << GetName() << endl;
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);
EStatus status;
if( (status = THaNonTrackingDetector::Init( date )) )
return fStatus=status;
for(Int_t ip=0;ip<fNLayers;ip++) {
if((status = fPlanes[ip]->Init( date ))) {
return fStatus=status;
}
}
// 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.
// 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;
fNLayers = 4; // Hardwire for now
BlockThick = new Double_t [fNLayers];
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 [fNLayers];
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 [fNLayers];
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*fNLayers];
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 );
cout << "THcShower::DefineVariables called " << GetName() << endl;
// 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);
Int_t nexthit = 0;
for(Int_t ip=0;ip<fNLayers;ip++) {
nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
}
/*
// 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)
////////////////////////////////////////////////////////////////////////////////