Skip to content

Feature/acastill deconvolution database #683

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nathanielerowe merged 8 commits into from
Jun 27, 2025
Merged

Feature/acastill deconvolution database #683

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
54e0e96
Adapt PMT deconvolution to work with calibration database
asanchezcastillo Mar 3, 2025
e44159f
Fix deconvolution includes
asanchezcastillo Mar 5, 2025
a82ac50
Merge branch 'develop' into feature/acastill_deconvolution_database
asanchezcastillo Mar 14, 2025
d804d0e
Merge branch 'develop' into feature/acastill_deconvolution_database
asanchezcastillo Mar 26, 2025
067256f
Fix include
asanchezcastillo Mar 26, 2025
a38c953
Merge branch 'develop' into feature/acastill_deconvolution_database
nathanielerowe May 27, 2025
c1cd4d5
Merge branch 'develop' into feature/acastill_deconvolution_database
nathanielerowe Jun 27, 2025
8ac8efe
Merge branch 'develop' into feature/acastill_deconvolution_database
nathanielerowe Jun 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
89 changes: 26 additions & 63 deletions sbndcode/OpDetReco/OpDeconvolution/Alg/OpDeconvolutionAlgWienerData_tool.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,8 @@

#include "sbndcode/OpDetReco/OpDeconvolution/Alg/OpDeconvolutionAlg.hh"

#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h"
#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h"

namespace opdet {
class OpDeconvolutionAlgWiener;
Expand All @@ -37,7 +39,7 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg {
public:
explicit OpDeconvolutionAlgWiener(fhicl::ParameterSet const& p);

~OpDeconvolutionAlgWiener() {}
~OpDeconvolutionAlgWiener() {delete fFilterTF1;}

// Required functions.
std::vector<raw::OpDetWaveform> RunDeconvolution(std::vector<raw::OpDetWaveform> const& wfHandle) override;
Expand Down Expand Up @@ -65,7 +67,6 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg {
double fPMTChargeToADC;
double fDecoWaveformPrecision;
short unsigned int fBaselineSample;
std::string fOpDetDataFile;
std::string fFilter;
std::string fElectronics;
bool fScaleHypoSignal;
Expand All @@ -80,6 +81,7 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg {
unsigned int NDecoWf;

TF1 *fFilterTF1;

std::vector<double> fSignalHypothesis;
std::vector<double> fNoiseHypothesis;

Expand All @@ -88,6 +90,8 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg {
short unsigned int fBaseSampleBins;
double fBaseVarCut;

sbndDB::PMTCalibrationDatabase const* fPMTCalibrationDatabaseService;

// Declare member functions
void ApplyExpoAvSmoothing(std::vector<double>& wf);
void ApplyUnAvSmoothing(std::vector<double>& wf);
Expand Down Expand Up @@ -125,7 +129,6 @@ opdet::OpDeconvolutionAlgWiener::OpDeconvolutionAlgWiener(fhicl::ParameterSet co
fPMTChargeToADC = p.get< double >("PMTChargeToADC");
fDecoWaveformPrecision = p.get< double >("DecoWaveformPrecision");
fBaselineSample = p.get< short unsigned int >("BaselineSample");
fOpDetDataFile = p.get< std::string >("OpDetDataFile");
fSkipChannelList = p.get< std::vector<int>>("SkipChannelList");
fFilter = p.get< std::string >("Filter");
fElectronics = p.get< std::string >("Electronics");
Expand All @@ -138,6 +141,8 @@ opdet::OpDeconvolutionAlgWiener::OpDeconvolutionAlgWiener(fhicl::ParameterSet co
fBaseVarCut = p.get< double >("BaseVarCut");
fFilterParams = p.get< std::vector<double> >("FilterParams");

fFilterTF1 = new TF1("FilterTemplate", fFilter.c_str());

fNormUnAvSmooth=1./(2*fUnAvNeighbours+1);
NDecoWf=0;
MaxBinsFFT=std::pow(2, fMaxFFTSizePow);
Expand All @@ -147,38 +152,10 @@ opdet::OpDeconvolutionAlgWiener::OpDeconvolutionAlgWiener(fhicl::ParameterSet co
if (fElectronics=="Daphne") fSamplingFreq=fDaphne_Freq/1000.;//in GHz
auto const* lar_prop = lar::providerFrom<detinfo::LArPropertiesService>();

//Load SER
std::string fname;
cet::search_path sp("FW_SEARCH_PATH");
sp.find_file(fOpDetDataFile, fname);
TFile* file = TFile::Open(fname.c_str(), "READ");
std::vector<std::vector<double>>* SinglePEVec_p;
std::vector<int>* fSinglePEChannels_p;
std::vector<double>* fPeakAmplitude_p;
std::vector<std::vector<double>> * fFilterParamVector_p;

file->GetObject("SERChannels", fSinglePEChannels_p);
file->GetObject("SinglePEVec", SinglePEVec_p);
file->GetObject("PeakAmplitude", fPeakAmplitude_p);
file->GetObject("FilterParams", fFilterParamVector_p);

if (fElectronics=="Daphne") file->GetObject("SinglePEVec_40ftCable_Daphne", SinglePEVec_p);
fSinglePEWaveVector = *SinglePEVec_p;
fSinglePEChannels = *fSinglePEChannels_p;
fPeakAmplitude = *fPeakAmplitude_p;
fFilterParamVector = *fFilterParamVector_p;

mf::LogInfo("OpDeconvolutionAlg")<<"Loaded SER from "<<fOpDetDataFile<<"... size="<<fSinglePEWave.size()<<std::endl;
file->Close();

if(fUseParamFilter && !fUseParamFilterInidividualChannel){
//If use param filter but not one for each channel, build it here
fFilterTF1 = new TF1("FilterTemplate", fFilter.c_str());
for(size_t k=0; k<fFilterParams.size(); k++)
fFilterTF1->SetParameter(k, fFilterParams[k]);
mf::LogInfo("OpDeconvolutionAlg")<<"Creating parametrized filter... TF1:"<<fFilter<<std::endl;
}
else if (!fUseParamFilter){
//Load PMT Calibration Database
fPMTCalibrationDatabaseService = lar::providerFrom<sbndDB::IPMTCalibrationDatabaseService const>();

if (!fUseParamFilter){
//Create light signal hypothesis for "on-fly" Wiener filter
fSignalHypothesis.resize(MaxBinsFFT, 0);
if(fFilter=="Wiener")
Expand All @@ -200,38 +177,26 @@ std::vector<raw::OpDetWaveform> opdet::OpDeconvolutionAlgWiener::RunDeconvolutio
{
int channelNumber = wf.ChannelNumber();
auto it = std::find(fSkipChannelList.begin(), fSkipChannelList.end(), channelNumber);
bool AnalyseChannel = false;
if (it == fSkipChannelList.end()) {
//If it's not try to find its SER in the file
for(size_t i=0; i<fSinglePEChannels.size(); i++)
fSinglePEWave = fPMTCalibrationDatabaseService->getSER(channelNumber);
double SPEAmplitude = fPMTCalibrationDatabaseService->getSPEAmplitude(channelNumber);
double SPEPeakValue = *std::max_element(fSinglePEWave.begin(), fSinglePEWave.end(), [](double a, double b) {return std::abs(a) < std::abs(b);});
double SinglePENormalization = std::abs(SPEAmplitude/SPEPeakValue);
std::transform(fSinglePEWave.begin(), fSinglePEWave.end(), fSinglePEWave.begin(), [SinglePENormalization](double val) {return val * SinglePENormalization;});
fSinglePEWave.resize(MaxBinsFFT, 0);
// If use channel dependent param filter
if(fUseParamFilter)
{
if(fSinglePEChannels[i]==channelNumber)
double GaussFilterPower = fPMTCalibrationDatabaseService->getGaussFilterPower(channelNumber);
double GaussFilterWC = fPMTCalibrationDatabaseService->getGaussFilterWC(channelNumber);
fFilterParamChannel = {GaussFilterWC, GaussFilterPower};
for(size_t k=0; k<fFilterParamChannel.size(); k++)
{
fSinglePEWave = fSinglePEWaveVector[i];
double SPEPeakValue = *std::max_element(fSinglePEWave.begin(), fSinglePEWave.end(), [](double a, double b) {return std::abs(a) < std::abs(b);});
double SinglePENormalization = std::abs(fPeakAmplitude[i]/SPEPeakValue);
std::transform(fSinglePEWave.begin(), fSinglePEWave.end(), fSinglePEWave.begin(), [SinglePENormalization](double val) {return val * SinglePENormalization;});
fSinglePEWave.resize(MaxBinsFFT, 0);
AnalyseChannel = true;
// If use channel dependent param filter
if(fUseParamFilterInidividualChannel)
{
fFilterParamChannel = fFilterParamVector[i];
fFilterTF1 = new TF1("FilterTemplate", fFilter.c_str());
for(size_t k=0; k<fFilterParamChannel.size(); k++)
{
fFilterTF1->SetParameter(k, fFilterParamChannel[k]);
}


mf::LogInfo("OpDeconvolutionAlg")<<"Creating parametrized filter... TF1:"<<fFilter << " for channel " << channelNumber <<std::endl;
}
break;
fFilterTF1->SetParameter(k, fFilterParamChannel[k]);
}
mf::LogInfo("OpDeconvolutionAlg")<<"Creating parametrized filter... TF1:"<<fFilter << " for channel " << channelNumber <<std::endl;
}
if(AnalyseChannel == false) mf::LogError("OpDeconvolutionAlg") << " SER for channel " << channelNumber <<" not found in the file \n";
}
if(!AnalyseChannel) continue;
//Read waveform
size_t wfsize=wf.Waveform().size();
if(wfsize>MaxBinsFFT){
Expand Down Expand Up @@ -537,15 +502,13 @@ std::vector<TComplex> opdet::OpDeconvolutionAlgWiener::DeconvolutionKernel(size_
}
}


if(fDebug){
std::string name="h_wienerfilter_"+std::to_string(NDecoWf);
TH1F * hs_wiener = tfs->make< TH1F >
(name.c_str(),"Wiener Filter;Frequency Bin;Magnitude",size/2, 0, size/2);
for(size_t k=0; k<size/2; k++)
hs_wiener->SetBinContent(k, TComplex::Abs( kernel[k]*serfft[k] ) );
}
if(fUseParamFilter && fUseParamFilterInidividualChannel) delete fFilterTF1;
return kernel;
}

Expand Down
1 change: 0 additions & 1 deletion sbndcode/OpDetReco/OpDeconvolution/Alg/opdeconvolution_alg_data.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@ OpDeconvolutionAlgData:
tool_type: "OpDeconvolutionAlgWienerData"
Debug: false
MaxFFTSizePow: 16
OpDetDataFile: "./OpDetSim/digi_pmt_sbnd_data_OV6.root"
BaseSampleBins: 30
BaseVarCut: 50
SkipChannelList: []
Expand Down