diff --git a/sbndcode/OpDetReco/OpDeconvolution/Alg/OpDeconvolutionAlgWienerData_tool.cc b/sbndcode/OpDetReco/OpDeconvolution/Alg/OpDeconvolutionAlgWienerData_tool.cc index cbd23acef..d33cc57ee 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/Alg/OpDeconvolutionAlgWienerData_tool.cc +++ b/sbndcode/OpDetReco/OpDeconvolution/Alg/OpDeconvolutionAlgWienerData_tool.cc @@ -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; @@ -37,7 +39,7 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg { public: explicit OpDeconvolutionAlgWiener(fhicl::ParameterSet const& p); - ~OpDeconvolutionAlgWiener() {} + ~OpDeconvolutionAlgWiener() {delete fFilterTF1;} // Required functions. std::vector RunDeconvolution(std::vector const& wfHandle) override; @@ -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; @@ -80,6 +81,7 @@ class opdet::OpDeconvolutionAlgWiener : opdet::OpDeconvolutionAlg { unsigned int NDecoWf; TF1 *fFilterTF1; + std::vector fSignalHypothesis; std::vector fNoiseHypothesis; @@ -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& wf); void ApplyUnAvSmoothing(std::vector& wf); @@ -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>("SkipChannelList"); fFilter = p.get< std::string >("Filter"); fElectronics = p.get< std::string >("Electronics"); @@ -138,6 +141,8 @@ opdet::OpDeconvolutionAlgWiener::OpDeconvolutionAlgWiener(fhicl::ParameterSet co fBaseVarCut = p.get< double >("BaseVarCut"); fFilterParams = p.get< std::vector >("FilterParams"); + fFilterTF1 = new TF1("FilterTemplate", fFilter.c_str()); + fNormUnAvSmooth=1./(2*fUnAvNeighbours+1); NDecoWf=0; MaxBinsFFT=std::pow(2, fMaxFFTSizePow); @@ -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(); - //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>* SinglePEVec_p; - std::vector* fSinglePEChannels_p; - std::vector* fPeakAmplitude_p; - std::vector> * 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 "<SetParameter(k, fFilterParams[k]); - mf::LogInfo("OpDeconvolutionAlg")<<"Creating parametrized filter... TF1:"<(); + + if (!fUseParamFilter){ //Create light signal hypothesis for "on-fly" Wiener filter fSignalHypothesis.resize(MaxBinsFFT, 0); if(fFilter=="Wiener") @@ -200,38 +177,26 @@ std::vector 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; igetSER(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; kSetParameter(k, fFilterParamChannel[k]); - } - - - mf::LogInfo("OpDeconvolutionAlg")<<"Creating parametrized filter... TF1:"<SetParameter(k, fFilterParamChannel[k]); } + mf::LogInfo("OpDeconvolutionAlg")<<"Creating parametrized filter... TF1:"<MaxBinsFFT){ @@ -537,7 +502,6 @@ std::vector opdet::OpDeconvolutionAlgWiener::DeconvolutionKernel(size_ } } - if(fDebug){ std::string name="h_wienerfilter_"+std::to_string(NDecoWf); TH1F * hs_wiener = tfs->make< TH1F > @@ -545,7 +509,6 @@ std::vector opdet::OpDeconvolutionAlgWiener::DeconvolutionKernel(size_ for(size_t k=0; kSetBinContent(k, TComplex::Abs( kernel[k]*serfft[k] ) ); } - if(fUseParamFilter && fUseParamFilterInidividualChannel) delete fFilterTF1; return kernel; } diff --git a/sbndcode/OpDetReco/OpDeconvolution/Alg/opdeconvolution_alg_data.fcl b/sbndcode/OpDetReco/OpDeconvolution/Alg/opdeconvolution_alg_data.fcl index fb1f68deb..6d6821780 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/Alg/opdeconvolution_alg_data.fcl +++ b/sbndcode/OpDetReco/OpDeconvolution/Alg/opdeconvolution_alg_data.fcl @@ -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: []