-
Vardan Tadevosyan authoredVardan Tadevosyan authored
THcShowerArray.cxx 11.16 KiB
//*-- Author :
//////////////////////////////////////////////////////////////////////////
//
// THcShowerArray
//
//////////////////////////////////////////////////////////////////////////
#include "THcShowerArray.h"
#include "TClonesArray.h"
#include "THcSignalHit.h"
#include "THcGlobals.h"
#include "THcParmList.h"
#include "THcHitList.h"
#include "THcShower.h"
#include "THcRawShowerHit.h"
#include "TClass.h"
#include "math.h"
#include "THaTrack.h"
#include "THaTrackProj.h"
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <fstream>
using namespace std;
ClassImp(THcShowerArray)
//______________________________________________________________________________
THcShowerArray::THcShowerArray( const char* name,
const char* description,
const Int_t layernum,
THaDetectorBase* parent )
: THaSubDetector(name,description,parent)
{
fADCHits = new TClonesArray("THcSignalHit",100);
fLayerNum = layernum;
}
THcShowerArray::~THcShowerArray()
{
// Destructor
delete fADCHits;
delete [] fA;
delete [] fP;
}
//_____________________________________________________________________________
THaAnalysisObject::EStatus THcShowerArray::Init( const TDatime& date )
{
// Extra initialization for shower layer: set up DataDest map
if( IsZombie())
return fStatus = kInitError;
// How to get information for parent
// if( GetParent() )
// fOrigin = GetParent()->GetOrigin();
EStatus status;
if( (status=THaSubDetector::Init( date )) )
return fStatus = status;
return fStatus = kOK;
}
//_____________________________________________________________________________
Int_t THcShowerArray::ReadDatabase( const TDatime& date )
{
char prefix[2];
prefix[0]=tolower(GetParent()->GetPrefix()[0]);
prefix[1]='\0';
cout << "Parent name: " << GetParent()->GetPrefix() << endl;
fUsingFADC=0;
fPedSampLow=0;
fPedSampHigh=9;
fDataSampLow=23;
fDataSampHigh=49;
DBRequest list[]={
{"cal_nrows", &fNRows, kInt},
{"cal_ncolumns", &fNColumns, kInt},
{"cal_using_fadc", &fUsingFADC, kInt, 0, 1},
{"cal_ped_sample_low", &fPedSampLow, kInt, 0, 1},
{"cal_ped_sample_high", &fPedSampHigh, kInt, 0, 1},
{"cal_data_sample_low", &fDataSampLow, kInt, 0, 1},
{"cal_data_sample_high", &fDataSampHigh, kInt, 0, 1},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list, prefix);
fNelem = fNRows*fNColumns;
// Here read the 2-D arrays of pedestals, gains, etc.
// Pedestal limits per channel.
fPedLimit = new Int_t [fNelem];
THcShower* fParent;
fParent = (THcShower*) GetParent();
for(Int_t i=0;i<fNelem;i++) {
fPedLimit[i] = fParent->GetPedLimit(i,fLayerNum,0); //layer 2, neg. side
}
fMinPeds = fParent->GetMinPeds();
InitializePedestals();
// Event by event amplitude and pedestal
fA = new Double_t[fNelem];
fP = new Double_t[fNelem];
#ifdef HITPIC
hitpic = new char*[fNRows];
for(Int_t row=0;row<fNRows;row++) {
hitpic[row] = new char[NPERLINE*(fNColumns+1)+2];
}
piccolumn=0;
#endif
// Debug output.
if (fParent->fdbg_init_cal) {
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerArray::ReadDatabase for "
<< GetParent()->GetPrefix() << ":" << endl;
cout << " Layer #" << fLayerNum << ", number of elements " << fNelem
<< endl;
cout << " Columns " << fNColumns << ", Rows " << fNRows << endl;
cout << " Using FADC " << fUsingFADC << endl; if (fUsingFADC) {
cout << " FADC pedestal sample low = " << fPedSampLow << ", high = "
<< fPedSampHigh << endl;
cout << " FADC data sample low = " << fDataSampLow << ", high = "
<< fDataSampHigh << endl;
}
// cout << " Origin of Layer at X = " << fOrigin.X()
// << " Y = " << fOrigin.Y() << " Z = " << fOrigin.Z() << endl;
cout << " fPedLimit:";
for(Int_t i=0;i<fNelem;i++) cout << " " << fPedLimit[i];
cout << endl;
cout << " fMinPeds = " << fMinPeds << endl;
cout << "---------------------------------------------------------------\n";
}
return kOK;
}
//_____________________________________________________________________________
Int_t THcShowerArray::DefineVariables( EMode mode )
{
// Initialize global variables
if( mode == kDefine && fIsSetup ) return kOK;
fIsSetup = ( mode == kDefine );
// Register variables in global list
RVarDef vars[] = {
{"adchits", "List of ADC hits", "fADCHits.THcSignalHit.GetPaddleNumber()"},
{"a", "Raw ADC Amplitude", "fA"},
{"p", "Dynamic ADC Pedestal", "fP"},
{ 0 }
};
return DefineVarsFromList( vars, mode );
}
//_____________________________________________________________________________
void THcShowerArray::Clear( Option_t* )
{
// Clears the hit lists
fADCHits->Clear();
}
//_____________________________________________________________________________
Int_t THcShowerArray::Decode( const THaEvData& evdata )
{
// Doesn't actually get called. Use Fill method instead
return 0;
}
//_____________________________________________________________________________
Int_t THcShowerArray::CoarseProcess( TClonesArray& tracks )
{
// Nothing is done here. See ProcessHits method instead.
//
return 0;
}
//_____________________________________________________________________________
Int_t THcShowerArray::FineProcess( TClonesArray& tracks )
{
return 0;}
//_____________________________________________________________________________
Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
{
// Extract the data for this layer from hit list
// Assumes that the hit list is sorted by layer, so we stop when the
// plane doesn't agree and return the index for the next hit.
THcShower* fParent;
fParent = (THcShower*) GetParent();
// Initialize variables.
fADCHits->Clear();
for(Int_t i=0;i<fNelem;i++) {
fA[i] = 0;
}
// Process raw hits. Get ADC hits for the plane, assign variables for each
// channel.
Int_t nrawhits = rawhits->GetLast()+1;
Int_t ihit = nexthit;
Int_t ngood = 0;
Int_t threshold = 100;
while(ihit < nrawhits) {
THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
// Should probably check that counter # is in range
if(fUsingFADC) {
fA[hit->fCounter-1] = hit->GetData(0,fPedSampLow,fPedSampHigh,
fDataSampLow,fDataSampHigh);
fP[hit->fCounter-1] = hit->GetPedestal(0,fPedSampLow,fPedSampHigh);
} else {
fA[hit->fCounter-1] = hit->GetData(0);
}
if(fA[hit->fCounter-1] > threshold) {
ngood++;
}
// Do other stuff like comparison to thresholds, signal hits, energy sums
ihit++;
}
#if 0
if(ngood > 0) {
cout << "+";
for(Int_t column=0;column<fNColumns;column++) {
cout << "-";
}
cout << "+" << endl;
for(Int_t row=0;row<fNRows;row++) {
cout << "|";
for(Int_t column=0;column<fNColumns;column++) {
Int_t counter = column*fNRows + row;
if(fA[counter]>threshold) {
cout << "X";
} else {
cout << " ";
}
}
cout << "|" << endl;
}
}
#endif#ifdef HITPIC
if(ngood > 0) {
for(Int_t row=0;row<fNRows;row++) {
if(piccolumn==0) {
hitpic[row][0] = '|';
}
for(Int_t column=0;column<fNColumns;column++) {
Int_t counter = column*fNRows+row;
if(fA[counter] > threshold) {
hitpic[row][piccolumn*(fNColumns+1)+column+1] = 'X';
} else {
hitpic[row][piccolumn*(fNColumns+1)+column+1] = ' ';
}
hitpic[row][(piccolumn+1)*(fNColumns+1)] = '|';
}
}
piccolumn++;
if(piccolumn==NPERLINE) {
cout << "+";
for(Int_t pc=0;pc<NPERLINE;pc++) {
for(Int_t column=0;column<fNColumns;column++) {
cout << "-";
}
cout << "+";
}
cout << endl;
for(Int_t row=0;row<fNRows;row++) {
hitpic[row][(piccolumn+1)*(fNColumns+1)+1] = '\0';
cout << hitpic[row] << endl;
}
piccolumn = 0;
}
}
#endif
return(ihit);
}
//_____________________________________________________________________________
Int_t THcShowerArray::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
{
// Extract data for this plane from hit list and accumulate in
// arrays for subsequent pedestal calculations.
Int_t nrawhits = rawhits->GetLast()+1;
Int_t ihit = nexthit;
while(ihit < nrawhits) {
THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
// OK for hit list sorted by layer.
// if(hit->fPlane > fLayerNum) {
// break;
// }
Int_t element = hit->fCounter - 1; // Should check if in range
Int_t adc = fUsingFADC ?
hit->GetData(0,fPedSampLow,fPedSampHigh,fDataSampLow,fDataSampHigh)
:
hit->GetData(0);
if(adc <= fPedLimit[element]) {
fPedSum[element] += adc;
fPedSum2[element] += adc*adc;
fPedCount[element]++;
if(fPedCount[element] == fMinPeds/5) {
fPedLimit[element] = 100 + fPedSum[element]/fPedCount[element]; }
}
ihit++;
}
fNPedestalEvents++;
// Debug output.
if ( ((THcShower*) GetParent())->fdbg_raw_cal ) {
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerArray::AcculatePedestals for "
<< GetParent()->GetPrefix() << ":" << endl;
cout << "Processed hit list for plane " << GetName() << ":\n";
for (Int_t ih=nexthit; ih<nrawhits; ih++) {
THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ih);
Int_t adc = fUsingFADC ?
hit->GetData(0,fPedSampLow,fPedSampHigh,fDataSampLow,fDataSampHigh)
:
hit->GetData(0);
cout << " hit " << ih << ":"
<< " plane = " << hit->fPlane
<< " counter = " << hit->fCounter
<< " ADC = " << adc
<< endl;
}
cout << "---------------------------------------------------------------\n";
}
return(ihit);
}
//_____________________________________________________________________________
void THcShowerArray::CalculatePedestals( )
{
// Use the accumulated pedestal data to calculate pedestals.
for(Int_t i=0; i<fNelem;i++) {
fPed[i] = ((Float_t) fPedSum[i]) / TMath::Max(1, fPedCount[i]);
fSig[i] = sqrt(((Float_t)fPedSum2[i])/TMath::Max(1, fPedCount[i])
- fPed[i]*fPed[i]);
fThresh[i] = fPed[i] + TMath::Min(50., TMath::Max(10., 3.*fSig[i]));
}
// Debug output.
if ( ((THcShower*) GetParent())->fdbg_raw_cal ) {
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShowerArray::CalculatePedestals for"
<< GetParent()->GetPrefix() << ":" << endl;
cout << " ADC pedestals and thresholds for calorimeter plane "
<< GetName() << endl;
for(Int_t i=0; i<fNelem;i++) {
cout << " element " << i << ": "
<< " Pedestal = " << fPed[i]
<< " /threshold = " << fThresh[i]
<< endl;
}
cout << "---------------------------------------------------------------\n";
}
}
//_____________________________________________________________________________
void THcShowerArray::InitializePedestals( )
{
fNPedestalEvents = 0;
fPedSum = new Int_t [fNelem];
fPedSum2 = new Int_t [fNelem];
fPedCount = new Int_t [fNelem];
fSig = new Float_t [fNelem];
fPed = new Float_t [fNelem];
fThresh = new Float_t [fNelem];
for(Int_t i=0;i<fNelem;i++) {
fPedSum[i] = 0;
fPedSum2[i] = 0;
fPedCount[i] = 0;
}
}