diff --git a/duneopdet/OpticalDetector/Deconvolution/DecoAnalysis_module.cc b/duneopdet/OpticalDetector/Deconvolution/DecoAnalysis_module.cc index af4af55f..d63b008e 100644 --- a/duneopdet/OpticalDetector/Deconvolution/DecoAnalysis_module.cc +++ b/duneopdet/OpticalDetector/Deconvolution/DecoAnalysis_module.cc @@ -3,12 +3,17 @@ // This analyzer module for Deconvolution creates histograms // with information from OpDetWaveforms and OpWaveforms. // @authors : Daniele Guffanti, Maritza Delgado, Sergio Manthey Corchado -// @created : 2022 -//=========================================================================== +// @created : 2022 +// +// Modified: Oct 7, 2024, Viktor Pec +// Added raw waveforms to the output. Redid directory structure +// inside the root file. Set the maximum number of histograms to +// save to 400. +// =========================================================================== #ifndef DecoAnalysis_h #define DecoAnalysis_h - + // LArSoft includes #include "larcore/Geometry/Geometry.h" #include "lardataobj/RawData/OpDetWaveform.h" @@ -24,11 +29,12 @@ #include "art/Framework/Principal/Handle.h" #include "canvas/Persistency/Common/Ptr.h" #include "canvas/Persistency/Common/PtrVector.h" +#include "canvas/Persistency/Common/FindOne.h" #include "art/Framework/Services/Registry/ServiceHandle.h" #include "art_root_io/TFileService.h" #include "art_root_io/TFileDirectory.h" #include "messagefacility/MessageLogger/MessageLogger.h" - + // ROOT includes #include "TH1.h" #include "THStack.h" @@ -36,7 +42,7 @@ #include "TLorentzVector.h" #include "TVector3.h" #include "TTree.h" - + // C++ Includes #include #include @@ -48,102 +54,136 @@ namespace opdet { class DecoAnalysis : public art::EDAnalyzer{ - public: - // Standard constructor and destructor for an ART module. - DecoAnalysis(const fhicl::ParameterSet&); - virtual ~DecoAnalysis(); - // The analyzer routine, called once per event. - void analyze (const art::Event&); - - private: - // The parameters we'll read from the .fcl file. - std::string fInputModuleDeco; // Input tag for OpDet collection - std::string fInputModuleDigi; // Input tag for OpDet collection - std::string fInstanceName; // Input tag for OpDet collection - double fSampleFreq; // in MHz + public: + // Standard constructor and destructor for an ART module. + DecoAnalysis(const fhicl::ParameterSet&); + virtual ~DecoAnalysis(); + // The analyzer routine, called once per event. + void analyze (const art::Event&); + + private: + // The parameters we'll read from the .fcl file. + std::string fInputModuleDeco; // Input tag for OpDet collection + std::string fInputModuleDigi; // Input tag for OpDet collection + std::string fInstanceNameDeco; // Input tag for OpDet collection + std::string fInstanceNameDigi; // Input tag for OpDet collection + double fSampleFreq; // in MHz //float fTimeBegin; // in us // float fTimeEnd; // in us - }; -} -#endif + static const size_t MAX_WFMS = 400; + }; +} + +#endif namespace opdet { DEFINE_ART_MODULE(DecoAnalysis) } - + namespace opdet { //----------------------------------------------------------------------- // Constructor DecoAnalysis::DecoAnalysis(fhicl::ParameterSet const& pset) - : EDAnalyzer(pset){ + : EDAnalyzer(pset){ // Read the fcl-file fInputModuleDeco = pset.get< std::string >("InputModuleDeco"); fInputModuleDigi = pset.get< std::string >("InputModuleDigi"); - fInstanceName = pset.get< std::string >("InstanceName"); - + fInstanceNameDeco = pset.get< std::string >("InstanceNameDeco"); + fInstanceNameDigi = pset.get< std::string >("InstanceNameDigi"); + // Obtain parameters from DetectorClocksService - auto const clockData = art::ServiceHandle()->DataForJob(); + auto const clockData = art::ServiceHandle()->DataForJob(); fSampleFreq = clockData.OpticalClock().Frequency(); art::ServiceHandle< art::TFileService > tfs; } - + //----------------------------------------------------------------------- // Destructor - DecoAnalysis::~DecoAnalysis() { + DecoAnalysis::~DecoAnalysis() { } - + //----------------------------------------------------------------------- - void DecoAnalysis::analyze(const art::Event& evt){ + void DecoAnalysis::analyze(const art::Event& evt){ + auto mfi = mf::LogInfo("DecoAnalysis::analyze()"); + + // Map to store how many waveforms are on one optical channel std::map< int, int > mapChannelWF; - + // Get deconvolved "OpWaveforms" from the event art::Handle< std::vector< recob::OpWaveform >> deconv; - evt.getByLabel(fInputModuleDeco, fInstanceName, deconv); + evt.getByLabel(fInputModuleDeco, fInstanceNameDeco, deconv); std::vector< recob::OpWaveform > OpWaveform; for (auto const& wf : *deconv) {OpWaveform.push_back(wf);} - + // Get Waveforms "OpDetWaveforms" from the event - art::Handle< std::vector< raw::OpDetWaveform >> digi; - evt.getByLabel(fInputModuleDigi, fInstanceName, digi); - std::vector< raw::OpDetWaveform > OpDetWaveform; - for (auto const& wf : *digi) {OpDetWaveform.push_back(wf);} + //art::Handle< std::vector< raw::OpDetWaveform >> digi; + art::FindOne fopDigi(deconv, evt, fInputModuleDeco); // Access ART's TFileService, which will handle creating and writing histograms for us art::ServiceHandle< art::TFileService > tfs; - + auto evt_dir = tfs->mkdir(TString::Format("run_%d_evt_%d", evt.id().run(), evt.id().event()).Data()); + double firstWaveformTime = -9999; - for (auto const& waveform : *digi) - { - if (firstWaveformTime < 0) firstWaveformTime = waveform.TimeStamp(); - } + for (long unsigned int i = 0; i < fopDigi.size(); ++i) + { + auto const &waveform = fopDigi.at(i).ref(); + if (firstWaveformTime == -9999 || firstWaveformTime > waveform.TimeStamp()) + firstWaveformTime = waveform.TimeStamp(); + } + + std::map > chan_dir_map; + + + mfi<<"\n" + <<"Number of deconvolved waveforms in this event: "< OpWaveform.size()) + max_wfms = OpWaveform.size(); + + mfi<<"Number of deconvolved waveforms to be saved in this event: "<make< TH1D >(histName.str().c_str(), TString::Format(";t - %f (#mus);",firstWaveformTime), decowaveform.NSignal(), startTime, endTime); + TH1D *decowaveformHist = chan_dir.second.make< TH1D >(TString::Format("waveform_%d", nth_waveform), TString::Format(";t - %f (#mus);",firstWaveformTime), decowaveform.NSignal(), startTime, endTime); // Copy values from the waveform into the histogram for(size_t tick = 0; tick < decowaveform.NSignal(); tick++){ - // Fill histogram with waveform - decowaveformHist->SetBinContent(tick, decowaveform.Signal()[tick]); - } - } - } + // Fill histogram with waveform + decowaveformHist->SetBinContent(tick+1, decowaveform.Signal()[tick]); + } + + // Create a new histogram with raw waveform + TH1D *rawwaveformHist = chan_dir.first.make< TH1D >(TString::Format("waveform_%d", nth_waveform), TString::Format(";t - %f (#mus);",firstWaveformTime), waveform.Waveform().size(), startTime, endTime); + // Copy values from the waveform into the histogram + for(size_t tick = 0; tick < waveform.Waveform().size(); tick++){ + // Fill histogram with waveform + rawwaveformHist->SetBinContent(tick+1, waveform.Waveform()[tick]); + } + } + } } // namespace opdet diff --git a/duneopdet/OpticalDetector/Deconvolution/Deconvolution.fcl b/duneopdet/OpticalDetector/Deconvolution/Deconvolution.fcl index 7047d6ea..e81dcc7c 100644 --- a/duneopdet/OpticalDetector/Deconvolution/Deconvolution.fcl +++ b/duneopdet/OpticalDetector/Deconvolution/Deconvolution.fcl @@ -1,3 +1,5 @@ +#include "dune_opdet_channels.fcl" + BEGIN_PROLOG ################################################################### @@ -8,48 +10,118 @@ WfmWienerfilter: { Name: "Wiener" Cutoff: 1. # Cutoff is not used in this filter } - + WfmGaussfilter: { Name: "Gauss" Cutoff: 0.13 # In MHz.The cutoff frequency is defined by the standard deviation in the frequency domain. } # The cutoff value should be changed if signal smoothing is not observed. - + ################################################################### -dune_deconvolution: +generic_dune_deconvolution: { module_type: "Deconvolution" - InputModule: "opdigi" - InstanceName: "" - + InputModule: "opdigi" + InstanceName: "" + #The LineNoiseRMS,PreTrigger,Pedestal, Samples and DigiDataFile below have the same values as the Digitizer. LineNoiseRMS: 2.6 # Pedestal RMS in [ADC] counts, likely an underestimate - PreTrigger: 100 # In [ticks] + PreTrigger: 100 # In [ticks] Pedestal: 1500 # In [ADC] Samples: 1000 # Timewindow (ReadoutWindow) in [ticks] Scale: 1 # Scaling of resulting deconvolution signal. - DigiDataColumn: 0 # SPE template source file column. - DigiDataFile: "SPE_DAPHNE2_FBK_2022.dat" - # The SPE template with undershoot and without pretrigger (in ADC*us), - + SPETemplateFileDataColumn: 0 # SPE template source file column. + SPETemplateFiles: [] # The SPE template with undershoot and without pretrigger (in ADC*us), + SPETemplateMapChannels: [] + SPETemplateMapTemplates: [] + IgnoreChannels: [] + + NoiseTemplateFiles: [] + NoiseTemplateMapChannels: [] + NoiseTemplateMapTemplates: [] + + InputPolarity: 1 # polarity of the input raw waveform + AutoScale: true # Scaling based on SPE amplitude from template (Use "true" for Wiener Filter and # "false" for Gauss Filter). If set to false the value of Scale is used. ApplyPostBLCorrection: false # Correct baseline after the deconvolution process(Use "false" for Wiener). PedestalBuffer: 20 # In [ticks], should always be smaller than PreTrigger. ApplyPostfilter: true # Filter the waveforms after "Wiener" deconvolution. - + WfmFilter: @local::WfmWienerfilter # Write the filter: "WfmWienerfilter" or "WfmGaussfilter" # If gauss filter is used, the following parameters should be changed: #WfmGaussfilter.Cutoff: 0.13; AutoScale: false; Scale: 1.; ApplyPostBLCorrection: true; ApplyPostfilter: false; - WfmPostfilter: @local::WfmGaussfilter # Only available "Gauss" postfilter. + WfmPostfilter: @local::WfmGaussfilter # Only available "Gauss" postfilter. } - -#Postfilter Cutoff + +#Postfilter Cutoff dune_deconvolution.WfmPostfilter.Cutoff: 1.8 # Use this value only for postfilter. + +# For backwards compatibility. This table is currently used in many workflows for different detectors, data and MC. +dune_deconvolution: @local::generic_dune_deconvolution +dune_deconvolution.SPETemplateFiles: ["SPE_DAPHNE2_FBK_2022.dat"] +# Set SPE template channel map: channel = -1 for all channels +dune_deconvolution.SPETemplateMapChannels: [-1] +dune_deconvolution.SPETemplateMapTemplates: [ 0] + + +##### Below, there are detector and data/MC spicific configurations. ##### + + +### FD ### +dune_fd_deconvolution: @local::dune_deconvolution + + #By debbuging and review the values (SNR,H,S,N,G0,G1,G,V,v) of the Wiener filter: -deconvolution_snr: @local::dune_deconvolution +deconvolution_snr: @local::dune_fd_deconvolution deconvolution_snr.OutputProduct: "SNR" -END_PROLOG \ No newline at end of file + +### ProtoDUNE HD ### +protodunehd_deconvolution: @local::dune_deconvolution +protodunehd_deconvolution.Samples: 1024 +protodunehd_deconvolution.Pedestal: 8180 +protodunehd_deconvolution.PreTrigger: 50 +protodunehd_deconvolution.PedestalBuffer: 30 +protodunehd_deconvolution.InputModule: "pdhddaphne" +protodunehd_deconvolution.InstanceName: "daq" +protodunehd_deconvolution.InputPolarity: -1 +protodunehd_deconvolution.ApplyPostfilter: true +protodunehd_deconvolution.WfmPostfilter.Cutoff: 1.5 +protodunehd_deconvolution.ApplyPostBLCorrection: true + +protodunehd_deconvolution.IgnoreChannels: @local::protodunehd_pds_channels.IgnoreChannels + +# SPE Template +protodunehd_deconvolution.SPETemplateFiles: @local::protodunehd_pds_channels.SPETemplateFiles +protodunehd_deconvolution.SPETemplateMapChannels: @local::protodunehd_pds_channels.SPETemplateMapChannels +protodunehd_deconvolution.SPETemplateMapTemplates: @local::protodunehd_pds_channels.SPETemplateMapTemplates +# Noise template +protodunehd_deconvolution.NoiseTemplateFiles: @local::protodunehd_pds_channels.NoiseTemplateFiles +protodunehd_deconvolution.NoiseTemplateMapChannels: @local::protodunehd_pds_channels.NoiseTemplateMapChannels +protodunehd_deconvolution.NoiseTemplateMapTemplates: @local::protodunehd_pds_channels.NoiseTemplateMapTemplates + + +### ProtoDUNE HD MC ### +protodunehd_mc_deconvolution: @local::protodunehd_deconvolution +# values below are as per settings for the digitizer,table protodunehd_opdigi from opticaldetectormodules_dune.fcl +protodunehd_mc_deconvolution.LineNoiseRMS: 4.5 +protodunehd_mc_deconvolution.Pedestal: 8200 +protodunehd_mc_deconvolution.PreTrigger: 128 +protodunehd_mc_deconvolution.InputModule: "opdigi" +protodunehd_mc_deconvolution.InstanceName: "" +protodunehd_mc_deconvolution.InputPolarity: 1 + +protodunehd_mc_deconvolution.IgnoreChannels: [] + +protodunehd_mc_deconvolution.SPETemplateFiles: @local::protodunehd_pds_channels_mc.SPETemplateFiles +protodunehd_mc_deconvolution.SPETemplateMapChannels: @local::protodunehd_pds_channels_mc.SPETemplateMapChannels +protodunehd_mc_deconvolution.SPETemplateMapTemplates: @local::protodunehd_pds_channels_mc.SPETemplateMapTemplates +# no noise templates used in PDHD MC +protodunehd_deconvolution.NoiseTemplateFiles: [] +protodunehd_deconvolution.NoiseTemplateMapChannels: [] +protodunehd_deconvolution.NoiseTemplateMapTemplates: [] + +END_PROLOG diff --git a/duneopdet/OpticalDetector/Deconvolution/Deconvolution_module.cc b/duneopdet/OpticalDetector/Deconvolution/Deconvolution_module.cc index b186eec9..ac1a1115 100644 --- a/duneopdet/OpticalDetector/Deconvolution/Deconvolution_module.cc +++ b/duneopdet/OpticalDetector/Deconvolution/Deconvolution_module.cc @@ -5,12 +5,16 @@ // Creating OpWaveform object // Filter wiener/gauss -FFT // @authors : Daniele Guffanti, Maritza Delgado, Sergio Manthey Corchado -// @created : Jan 26, 2022 -//============================================================================= - +// @created : Jan 26, 2022 +// +// Modified: +// Oct 7, 2024, Viktor Pec +// Added possibility to use SPE and noise templates by channel. +//============================================================================= + #ifndef Deconvolution_h #define Deconvolution_h - + // Framework includes #include "art/Framework/Core/EDProducer.h" @@ -30,11 +34,13 @@ #include "art/Framework/Services/Registry/ActivityRegistry.h" #include "messagefacility/MessageLogger/MessageLogger.h" +#include "lardata/Utilities/AssociationUtil.h" + // ART extensions #include "nurandom/RandomUtils/NuRandomService.h" // LArSoft includes -#include "larcore/CoreUtils/ServiceUtil.h" +#include "larcore/CoreUtils/ServiceUtil.h" #include "larcore/Geometry/Geometry.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" @@ -70,7 +76,7 @@ #include "TF1.h" #include "TComplex.h" #include "TFFTComplexReal.h" -#include +#include using std::string; using std::vector; @@ -81,32 +87,46 @@ namespace opdet { public: struct Config { public: - fhicl::Atom module_type{ fhicl::Name("module_type") }; - fhicl::Atom InputModule{ fhicl::Name("InputModule") }; - fhicl::Atom InstanceName{ fhicl::Name("InstanceName") }; + fhicl::Atom module_type{ fhicl::Name("module_type") }; + fhicl::Atom InputModule{ fhicl::Name("InputModule") }; + fhicl::Atom InstanceName{ fhicl::Name("InstanceName") }; fhicl::Atom LineNoiseRMS{ fhicl::Name("LineNoiseRMS"), 1.0 }; fhicl::Atom PreTrigger{ fhicl::Name("PreTrigger"), 0}; - fhicl::Atom Pedestal{ fhicl::Name("Pedestal"), 1500}; - fhicl::Atom DigiDataFile{ fhicl::Name("DigiDataFile") }; - fhicl::Atom DigiDataColumn{ fhicl::Name("DigiDataColumn"), 1 }; + fhicl::Atom Pedestal{ fhicl::Name("Pedestal"), 1500}; + fhicl::Sequence SPETemplateFiles{ fhicl::Name("SPETemplateFiles") }; + fhicl::Atom SPETemplateFileDataColumn{ fhicl::Name("SPETemplateFileDataColumn"), 1 }; + fhicl::Sequence NoiseTemplateFiles{ fhicl::Name("NoiseTemplateFiles") }; + fhicl::Atom Samples{ fhicl::Name("Samples"), 1000 }; - fhicl::Atom PedestalBuffer{ fhicl::Name("PedestalBuffer"), 10 }; - fhicl::Atom Scale{ fhicl::Name("Scale"), 1 }; + fhicl::Atom PedestalBuffer{ fhicl::Name("PedestalBuffer"), 10 }; + fhicl::Atom Scale{ fhicl::Name("Scale"), 1 }; fhicl::Atom ApplyPostfilter{ fhicl::Name("ApplyPostfilter"), false }; - fhicl::Atom ApplyPostBLCorrection{ fhicl::Name("ApplyPostBLCorrection") }; + fhicl::Atom ApplyPostBLCorrection{ fhicl::Name("ApplyPostBLCorrection") }; fhicl::Atom AutoScale{ fhicl::Name("AutoScale"), false }; - fhicl::Atom OutputProduct{ fhicl::Name("OutputProduct"), "decowave"}; + + fhicl::Atom InputPolarity{ fhicl::Name("InputPolarity") }; + + + fhicl::Atom OutputProduct{ fhicl::Name("OutputProduct"), "decowave"}; + + fhicl::Sequence SPETemplateMap_channels{ fhicl::Name("SPETemplateMapChannels") }; + fhicl::Sequence SPETemplateMap_templates{ fhicl::Name("SPETemplateMapTemplates") }; + fhicl::Sequence NoiseTemplateMap_channels{ fhicl::Name("NoiseTemplateMapChannels") }; + fhicl::Sequence NoiseTemplateMap_templates{ fhicl::Name("NoiseTemplateMapTemplates") }; + + fhicl::Sequence IgnoreChannels{ fhicl::Name("IgnoreChannels") }; // integer to allow for channel -1 = unrecognized channel + struct Filter { fhicl::Atom Name{fhicl::Name("Name"), "Gauss"}; - fhicl::Atom Cutoff{fhicl::Name("Cutoff"), 2}; + fhicl::Atom Cutoff{fhicl::Name("Cutoff"), 2}; }; fhicl::Table Filter{ fhicl::Name("WfmFilter") }; struct ExtraFilter { - fhicl::Atom Name{fhicl::Name("Name"), "Gauss"}; - fhicl::Atom Cutoff{fhicl::Name("Cutoff"), 2}; - }; + fhicl::Atom Name{fhicl::Name("Name"), "Gauss"}; + fhicl::Atom Cutoff{fhicl::Name("Cutoff"), 2}; + }; fhicl::Table Postfilter{ fhicl::Name("WfmPostfilter") }; }; @@ -115,11 +135,14 @@ namespace opdet { //! Input signal shape model,for Wiener filter could use delta or scintillation light signal. enum EInputShape {kDelta = 0, kScint = 1}; //! Waveform filter type - enum EFilterType {kOther = 0, kWiener = 1, kGauss = 2, kNone = 3}; + enum EFilterType {kOther = 0, kWiener = 1, kGauss = 2, kNone = 3}; + + enum {kDefaultChannel = -1}; + struct WfmExtraFilter_t { - TString fName; - Double_t fCutoff; + TString fName; + Double_t fCutoff; WfmExtraFilter_t() {fName = ""; fCutoff = 1.0;} @@ -129,86 +152,86 @@ namespace opdet { WfmExtraFilter_t(const struct Config::ExtraFilter& config) { fName = config.Name(); - fCutoff = config.Cutoff(); + fCutoff = config.Cutoff(); } - }; + }; struct WfmFilter_t { - TString fName; - EFilterType fType; - Double_t fCutoff; + TString fName; + EFilterType fType; + Double_t fCutoff; WfmFilter_t() : fName("Unknown_filter"), fType(kOther), fCutoff(32) {} WfmFilter_t(const char* name) : fName(name), fType(kOther), fCutoff(32) {} WfmFilter_t(const struct Config::Filter& config) { - fName = config.Name(); - fCutoff = config.Cutoff(); - if (fName == "Wiener") fType = kWiener; - else if (fName == "Gauss") fType = kGauss; - else fType = kOther; + fName = config.Name(); + fCutoff = config.Cutoff(); + if (fName == "Wiener") fType = kWiener; + else if (fName == "Gauss") fType = kGauss; + else fType = kOther; } }; /** * @brief Helper struct to handle waveforms' FT * - * Helper data struct to store and handle waveforms' Fourier Transform. + * Helper data struct to store and handle waveforms' Fourier Transform. */ struct CmplxWaveform_t { //! vector of Fourier coefficients - std::vector fCmplx; + std::vector fCmplx; //! vector of the real part of Fourier coefficients - std::vector fRe; + std::vector fRe; //! vector of the imaginary part of Fourier coefficients - std::vector fIm; + std::vector fIm; //! Empty constructor inline CmplxWaveform_t() {} //! Constructor inline CmplxWaveform_t(int size) { - fCmplx = std::vector(size, TComplex(0., 0)); - fRe = std::vector(size, 0.); - fIm = std::vector(size, 0.); + fCmplx = std::vector(size, TComplex(0., 0)); + fRe = std::vector(size, 0.); + fIm = std::vector(size, 0.); } //! Copy constructor inline CmplxWaveform_t(const CmplxWaveform_t& cwf) { - fCmplx = std::vector( cwf.fCmplx ); - fRe = std::vector( cwf.fRe ); - fIm = std::vector( cwf.fIm ); + fCmplx = std::vector( cwf.fCmplx ); + fRe = std::vector( cwf.fRe ); + fIm = std::vector( cwf.fIm ); } //! Destructor inline ~CmplxWaveform_t() { - fCmplx.clear(); fRe.clear(); fIm.clear(); + fCmplx.clear(); fRe.clear(); fIm.clear(); } /** * @brief Point-side product of complex waveforms * - * Compute the product of two complex waveforms point + * Compute the product of two complex waveforms point * by point */ inline CmplxWaveform_t operator*(const CmplxWaveform_t& cwf) const { - CmplxWaveform_t result(fCmplx.size()); + CmplxWaveform_t result(fCmplx.size()); if (cwf.fCmplx.size() != fCmplx.size()) { - printf("opdet::Deconvolution_module::CmplxWaveform_t::operator* ERROR"); - printf(" waveforms size does not match.\n"); + printf("opdet::Deconvolution_module::CmplxWaveform_t::operator* ERROR"); + printf(" waveforms size does not match.\n"); return result; } for (size_t i=0; i fSPETemplateFiles; //!< single p.e. template source file + size_t fSPETemplateFileDataColumn; //!< single p.e. template source file column + std::vector fNoiseTemplateFiles; //!< noise template source file + double fScale; //!< Scaling of resulting wvfs int fSamples; //!< Same value as ReadoutWindow in digitizer int fPedestalBuffer; //!< Used to calculate pedestal which is definded as PreTrigger - PedestalBuffer in [ticks] //!< default is 10 ticks -> should be adapted to resulting peak width (after deconvolution). - bool fApplyPostfilter; + bool fApplyPostfilter; bool fApplyPostBLCorr; bool fAutoScale; - + + short int fInputPolarity; //!< whether the input raw waveform is positive or negative + + + // Additional parameters here - + //--- Single photoelectron variables - std::vector fSinglePEWaveform; //!< Vector that stores the single PE template in ADC*us - double fSinglePEAmplitude; //!< single PE amplitude for found maximum peak in the template. + std::vector > fSinglePEWaveforms; //!< Vector that stores the single PE template in ADC*us + std::vector fSinglePEWaveforms_fft; //!< Fourier transform of the tamplates + std::vector fSinglePEAmplitudes; //!< single PE amplitude for found maximum peak in the template. unsigned int WfDeco; //!< Number of waveform processed - + std::map fChannelToTemplateMap; //!< maps a channel id to the input SPE template file (index in fSinglePEWaveforms) + unsigned short fUseSingleSPETemplate; + std::set fIgnoreChannels; //!< List of channels to ignore in deconvolution + + // Noise templates -- input in frequency domain + std::vector > fNoiseTemplates; //!< Vector that stores noise template in frequency domain + std::map fChannelToNoiseTemplateMap; //!< maps a channel id to the input SPE template file (index in fSingle + std::vector fNoiseDefault; + + //--------Filter Variables - std::string fOutputProduct; + std::string fOutputProduct; WfmExtraFilter_t fPostfilterConfig; - WfmFilter_t fFilterConfig; - + WfmFilter_t fFilterConfig; + //Input signal shape model for Wiener filter could use delta "kDelta" or scintillation light signals "kScint" - EInputShape fInputShape = kDelta; + EInputShape fInputShape = kDelta; // EInputShape fInputShape = kScint; //uncomment if set to kScint private: int CountFileColumns(const char* file_path); - void SourceSPEDigiDataFile(); - void BuildExtraFilter(CmplxWaveform_t& xF0, const WfmExtraFilter_t config); + void SourceSPETemplateFiles(); + void SourceNoiseTemplateFiles(); + void BuildExtraFilter(CmplxWaveform_t& xF0, const WfmExtraFilter_t config); void ComputeExpectedInput(std::vector& s, double nmax); - void CopyToOutput(const std::vector& v, std::vector& target); - void CopyToOutput(const CmplxWaveform_t& v, std::vector& target); - Double_t ComputeAutoNormalization(CmplxWaveform_t& xGH, const float thrs=0); - TVirtualFFT* fft_r2c; - TVirtualFFT* fft_c2r; - CmplxWaveform_t fxG0; - CmplxWaveform_t fxG1; + void CopyToOutput(const std::vector& v, std::vector& target); + void CopyToOutput(const std::vector& v, std::vector& target); + void CopyToOutput(const CmplxWaveform_t& v, std::vector& target); + Double_t ComputeAutoNormalization(CmplxWaveform_t& xGH, const float thrs=0); + TVirtualFFT* fft_r2c; + TVirtualFFT* fft_c2r; + CmplxWaveform_t fxG0; + CmplxWaveform_t fxG1; }; } #endif @@ -296,37 +336,46 @@ namespace opdet{ DEFINE_ART_MODULE(Deconvolution) } +namespace { + using AssnsRawDeconv = art::Assns< raw::OpDetWaveform, recob::OpWaveform >; +} + namespace opdet { //--------------------------------------------------------------------------- // Constructor Deconvolution::Deconvolution(const Parameters& pars) - : EDProducer{pars}, - fInputModule{ pars().InputModule()}, - fInstanceName{ pars().InstanceName()}, - - fPedestal{ pars().Pedestal()}, - fLineNoiseRMS{ pars().LineNoiseRMS() }, - fPreTrigger{ pars().PreTrigger()}, - - fDigiDataFile{ pars().DigiDataFile()}, - fDigiDataColumn{ pars().DigiDataColumn()}, - fScale{ pars().Scale()}, - fSamples{ pars().Samples()}, - fPedestalBuffer{ pars().PedestalBuffer()}, - fApplyPostfilter{ pars().ApplyPostfilter()}, - fApplyPostBLCorr{ pars().ApplyPostBLCorrection()}, - fAutoScale{ pars().AutoScale()}, - fOutputProduct{ pars().OutputProduct() }, - fPostfilterConfig{ WfmExtraFilter_t( pars().Postfilter()) }, - fFilterConfig{ WfmFilter_t( pars().Filter() ) }, - fxG0(fSamples), - fxG1(fSamples) - { - // Declare that we'll produce a vector of OpDetWaveforms - WfDeco=0; + : EDProducer{pars}, + fInputModule{ pars().InputModule()}, + fInstanceName{ pars().InstanceName()}, + fPedestal{ pars().Pedestal()}, + fLineNoiseRMS{ pars().LineNoiseRMS() }, + fPreTrigger{ pars().PreTrigger()}, + fSPETemplateFiles{ pars().SPETemplateFiles()}, + fSPETemplateFileDataColumn{ pars().SPETemplateFileDataColumn()}, + fNoiseTemplateFiles{ pars().NoiseTemplateFiles()}, + fScale{ pars().Scale()}, + fSamples{ pars().Samples()}, + fPedestalBuffer{ pars().PedestalBuffer()}, + fApplyPostfilter{ pars().ApplyPostfilter()}, + fApplyPostBLCorr{ pars().ApplyPostBLCorrection()}, + fAutoScale{ pars().AutoScale()}, + fInputPolarity{ pars().InputPolarity()}, + fUseSingleSPETemplate(0), + fNoiseDefault(fSamples, fLineNoiseRMS*fLineNoiseRMS*fSamples), + fOutputProduct{ pars().OutputProduct() }, + fPostfilterConfig{ WfmExtraFilter_t( pars().Postfilter()) }, + fFilterConfig{ WfmFilter_t( pars().Filter() ) }, + fxG0(fSamples), + fxG1(fSamples) + { + auto mfi = mf::LogInfo("Deconvolution::Deconvolution()"); + + // Declare that we'll produce a vector of OpDetWaveforms + WfDeco=0; // This module produces - produces< std::vector< recob::OpWaveform> > (); + produces< std::vector< recob::OpWaveform> > (); + produces< AssnsRawDeconv > (); // Obtain parameters from TimeService @@ -335,128 +384,240 @@ namespace opdet { fft_r2c = TVirtualFFT::FFT(1, &fSamples, "M R2C K"); fft_c2r = TVirtualFFT::FFT(1, &fSamples, "M C2R K"); - - SourceSPEDigiDataFile(); + + // Prepare the SPE waveform templates + SourceSPETemplateFiles(); + // Prepare the Fourier transforms + for (auto &xh: fSinglePEWaveforms) { + fSinglePEWaveforms_fft.push_back(CmplxWaveform_t(fSamples)); + auto &xH = fSinglePEWaveforms_fft.back(); + + fft_r2c->SetPoints(&xh[0]); + fft_r2c->Transform(); + fft_r2c->GetPointsComplex(&xH.fRe[0], &xH.fIm[0]); + xH.MakeCmplx(); + } + + + // prepare channel to template map + { + auto channels = pars().SPETemplateMap_channels(); + auto templates = pars().SPETemplateMap_templates(); + auto chann = channels.begin(); + auto templ = templates.begin(); + for (;chann != channels.end(); ++chann, ++templ) { + fChannelToTemplateMap[*chann] = *templ; + } + // deal with a case where a single channel-to-template map was given for all channels + if ( fChannelToTemplateMap.size() == 1 + && channels[0] == kDefaultChannel + && fSinglePEWaveforms.size() == 1) { + fUseSingleSPETemplate = 1; + } else { + // FIXME: Throw exception in case malformed configuration + } + + } + + // Prepare the noise templates + SourceNoiseTemplateFiles(); + { + auto channels = pars().NoiseTemplateMap_channels(); + auto templates = pars().NoiseTemplateMap_templates(); + auto chann = channels.begin(); + auto templ = templates.begin(); + for (;chann != channels.end(); ++chann, ++templ) { + fChannelToNoiseTemplateMap[*chann] = *templ; + } + } + + for ( auto chan : pars().IgnoreChannels() ) { + fIgnoreChannels.insert(chan); + } + + //=== info print out === + mfi<<"Input waveform polarity set to: " << fInputPolarity << "\n"; + // info on channel to SPE template map + mfi<<"Channels mapped to SPE template files:\n"; + { + std::map< std::string, std::vector > templ_to_channel_map; + for (auto itm: fChannelToTemplateMap) + templ_to_channel_map[fSPETemplateFiles[itm.second]].push_back(itm.first); + for (auto itm: templ_to_channel_map) { + mfi<<" "< > templ_to_channel_map; + for (auto itm: fChannelToNoiseTemplateMap) + templ_to_channel_map[fNoiseTemplateFiles[itm.second]].push_back(itm.first); + for (auto itm: templ_to_channel_map) { + mfi<<" "< > wfHandle; evt.getByLabel(fInputModule, fInstanceName, wfHandle); - + //****************************** //-- Read Waveform---- - //****************************** + //****************************** - std::vector digi_wave = *wfHandle; + std::vector digi_wave = *wfHandle; int NOpDetWaveform = digi_wave.size(); //Number of waveforms in OpDetWaveform Object - std::cout << NOpDetWaveform << std::endl; + mfi << "Number of waveforms to process: " << NOpDetWaveform << "\n"; //Pointer that will store produced DecoWaveform auto out_recob = std::make_unique< std::vector< recob::OpWaveform > >(); + // associations to the raw waveform + std::unique_ptr< AssnsRawDeconv > assnPtr(new AssnsRawDeconv); std::vector out_digiwave(fSamples); //vector in which the waveform will be saved std::vector out_recob_float(fSamples); //vector in which the decowaveform will be saved, using float - std::vector xv(fSamples, 0.); - - //****************************** - //--- Setup filter's components - //****************************** + std::vector xv(fSamples, 0.); + - // xH: FFT of spe response - std::vector xh(fSinglePEWaveform); - CmplxWaveform_t xH(fSamples); - fft_r2c->SetPoints(&xh[0]); - fft_r2c->Transform(); - fft_r2c->GetPointsComplex(&xH.fRe[0], &xH.fIm[0]); - xH.MakeCmplx(); - //****************************** //--- Process waveforms //****************************** + int iDigi = -1; // index in digital waveform vector for (auto const& wf: digi_wave) { - CmplxWaveform_t xV(fSamples); - CmplxWaveform_t xS(fSamples); - CmplxWaveform_t xN(fSamples); - CmplxWaveform_t xG(fSamples); - CmplxWaveform_t xY(fSamples); - CmplxWaveform_t xGH(fSamples); - std::vector xSNR(fSamples, 0.); + ++iDigi; + auto channel = wf.ChannelNumber(); + + // check if this channel is to be ignored + if ( fIgnoreChannels.count(channel) ) + continue; + + //auto &xh = fSinglePEWaveforms[fChannelToTemplateMap[channel]]; + // vvv allow to use default template for all channels vvv + auto effChannel = channel; + if (fUseSingleSPETemplate) + effChannel = kDefaultChannel; + // ^^^ + auto &xH = fSinglePEWaveforms_fft[fChannelToTemplateMap[effChannel]]; // get the SPE template relevant for this channel + auto &speapmlitude = fSinglePEAmplitudes[fChannelToTemplateMap[effChannel]]; + + CmplxWaveform_t xV(fSamples); + CmplxWaveform_t xS(fSamples); + CmplxWaveform_t xG(fSamples); + CmplxWaveform_t xY(fSamples); + CmplxWaveform_t xGH(fSamples); + std::vector xSNR(fSamples, 0.); int OriginalWaveformSize = wf.Waveform().size(); + + // noise + auto iTemplate = fChannelToNoiseTemplateMap.find(channel); + auto &xN = (iTemplate != fChannelToNoiseTemplateMap.end()) ? fNoiseTemplates[iTemplate->second] : fNoiseDefault; // get the noise template relevant for this channel + + //----------------------Resize deconvoluted signals (using floats) to original waveform size - if (static_cast(OriginalWaveformSize) <= fSamples) { - out_recob_float.resize(fSamples,0); + if (static_cast(OriginalWaveformSize) <= fSamples) { + out_recob_float.resize(fSamples,0); } else { out_recob_float.resize(fSamples); } + // Calculate pedestal + double pedestal = 0.; + for (int i = 0; i < int(fPreTrigger-fPedestalBuffer); ++i) + pedestal += wf[i]; + pedestal /= (fPreTrigger-fPedestalBuffer); + for (Int_t i= 0; i < fSamples; i++){ - // Remove baseline - if (i < static_cast(wf.Waveform().size())) xv[i] = (wf[i]-fPedestal); + // Remove baseline and deal with input waveform polarity: make sure xv has positive polarity + if (i < static_cast(wf.Waveform().size())) xv[i] = fInputPolarity*(wf[i]-pedestal); // if waveform is shorter than fSamples fill the rest with noise else xv[i] = CLHEP::RandGauss::shoot(0, fLineNoiseRMS); } - + //---------------------------------------------------- Guess input signal // Found maximum peak in the Waveform + // Assume xv has positive polarity Double_t SPE_Max = 0; double maxADC=*max_element(xv.begin(),xv.end()); double maxAmplit= maxADC; // Pedestal already subtracted - SPE_Max = maxAmplit/fSinglePEAmplitude; - + SPE_Max = maxAmplit/speapmlitude; + std::vectorxs(fSamples,0.); // Compute expected input (using a delta or the scint tile profile as a model) - ComputeExpectedInput(xs, SPE_Max); - + ComputeExpectedInput(xs, SPE_Max); + //-------------------------------------------------- Compute waveform FFT fft_r2c->SetPoints(&xv[0]); fft_r2c->Transform(); fft_r2c->GetPointsComplex(&xV.fRe[0], &xV.fIm[0]); - xV.MakeCmplx(); + xV.MakeCmplx(); //----------------------------------------------------- Compute input FFT fft_r2c->SetPoints(&xs[0]); fft_r2c->Transform(); fft_r2c->GetPointsComplex(&xS.fRe[0], &xS.fIm[0]); - xS.MakeCmplx(); + xS.MakeCmplx(); //****************************** // Compute filters. - //****************************** + //****************************** + + Double_t fFrequencyCutOff = fFilterConfig.fCutoff; + Double_t fTickCutOff = fSamples*fFrequencyCutOff/fSampleFreq; - Double_t fFrequencyCutOff = fFilterConfig.fCutoff; - Double_t fTickCutOff = fSamples*fFrequencyCutOff/fSampleFreq; - for (int i=0; iSetPointsComplex(&xY.fRe[0], &xY.fIm[0]); fft_c2r->Transform(); double *xy = fft_c2r->GetPointsReal(); - std::vector xvdec(xy, xy+fSamples); - + std::vector xvdec(xy, xy+fSamples); + Double_t scale = 1.0 / fSamples; - if (!fAutoScale) {scale = fScale / fSamples;} + if (!fAutoScale) {scale = fScale / fSamples;} else { // Apply filter to the detector response (for normalization) - xGH = xG * xH; - xGH.MakeReAndIm(); - Double_t filter_norm = ComputeAutoNormalization(xGH); + xGH = xG * xH; + xGH.MakeReAndIm(); + Double_t filter_norm = ComputeAutoNormalization(xGH); scale = filter_norm / (Double_t)fSamples; } - if (fApplyPostfilter) { - CmplxWaveform_t xxY(fSamples); - std::vector ytmp(xvdec.begin(), xvdec.end()); - fft_r2c->SetPoints(&ytmp[0]); - fft_r2c->Transform(); - fft_r2c->GetPointsComplex(&xxY.fRe[0], &xxY.fIm[0]); - xxY.MakeCmplx(); - xxY = xxY * fxG1; - xxY.MakeReAndIm(); - fft_c2r->SetPointsComplex(&xxY.fRe[0], &xxY.fIm[0]); - fft_c2r->Transform(); - double* xxy = fft_c2r->GetPointsReal(); - double g1_scale = 1.0 / fSamples; - for (int i=0; i ytmp(xvdec.begin(), xvdec.end()); + fft_r2c->SetPoints(&ytmp[0]); + fft_r2c->Transform(); + fft_r2c->GetPointsComplex(&xxY.fRe[0], &xxY.fIm[0]); + xxY.MakeCmplx(); + xxY = xxY * fxG1; + xxY.MakeReAndIm(); + fft_c2r->SetPointsComplex(&xxY.fRe[0], &xxY.fIm[0]); + fft_c2r->Transform(); + double* xxy = fft_c2r->GetPointsReal(); + double g1_scale = 1.0 / fSamples; + for (int i=0; iSetPointsComplex(&xV.fRe[0], &xV.fIm[0]); - fft_c2r->Transform(); - Double_t* vv = fft_c2r->GetPointsReal(); - for (int i=0; iSetPointsComplex(&xV.fRe[0], &xV.fIm[0]); + fft_c2r->Transform(); + Double_t* vv = fft_c2r->GetPointsReal(); + for (int i=0; iemplace_back(std::move(decwav)); + + out_recob->emplace_back(std::move(decwav)); + + // create association + art::Ptr opwfm_ptr(wfHandle, iDigi); + util::CreateAssn( *this, evt, *out_recob, opwfm_ptr, *(assnPtr.get()), out_recob->size()-1); }//waveforms loop - + //-------------------------------------Print fType Filter if (fFilterConfig.fType == Deconvolution::kWiener){printf("***Wiener Filter****");} if (fFilterConfig.fType == Deconvolution::kGauss){printf("***Gauss Filter***");} if (fFilterConfig.fType == Deconvolution::kOther){printf("***Standart dec***");} if (fApplyPostfilter){printf("***ApplyPostfilter***");} - + //------------------------------------------------ - + // Push the OpDetWaveforms and OpWaveform into the event evt.put(std::move(out_recob)); + evt.put(std::move(assnPtr)); WfDeco++; } @@ -578,50 +748,50 @@ namespace opdet { /** * @brief Build a filter to be applied after the deconvolution * - * Construct an extra filter to be applied after + * Construct an extra filter to be applied after * the deconvolution process. - * Different filters can be implemented by switching the flag + * Different filters can be implemented by switching the flag * `fPostfilterConfig.fName` - * via the Config::ExtraFilter::name parameter. + * via the Config::ExtraFilter::name parameter. * * @param xF */ void Deconvolution::BuildExtraFilter(CmplxWaveform_t& xF, const WfmExtraFilter_t config) { if (config.fName != "Gauss") { - printf("Deconvolution::BuildExtraFilter WARNING: Unknown filter model %s. Skip.\n", - config.fName.Data()); + printf("Deconvolution::BuildExtraFilter WARNING: Unknown filter model %s. Skip.\n", + config.fName.Data()); return; - } + } // Compute sigma corresponding to the given cutoff frequency - const Double_t df = fSampleFreq / (Double_t)fSamples; - const Double_t cutoff = config.fCutoff / df; - const Double_t k_cutoff = sqrt(log(2)); + const Double_t df = fSampleFreq / (Double_t)fSamples; + const Double_t cutoff = config.fCutoff / df; + const Double_t k_cutoff = sqrt(log(2)); const double sigma = fSamples * k_cutoff / (TMath::TwoPi() * cutoff); - const int mu = 0.5*fSamples; + const int mu = 0.5*fSamples; - printf("Deconvolution::BuildExtraFilter sigma is %g\n", sigma); + printf("Deconvolution::BuildExtraFilter sigma is %g\n", sigma); std::vector xf(fSamples, 0.); for (int i=0; i re_(fSamples, 0.); + std::vector re_(fSamples, 0.); std::vector im_(fSamples, 0.); - fft_r2c->SetPoints(&xf[0]); - fft_r2c->Transform(); + fft_r2c->SetPoints(&xf[0]); + fft_r2c->Transform(); fft_r2c->GetPointsComplex(&re_[0], &im_[0]); for (int i=0; i<0.5*fSamples+1; i++) { - TComplex F(re_.at(i), im_.at(i)); - TComplex phase = TComplex(0., -TMath::TwoPi()*i*mu/(fSamples)); - xF.fCmplx.at(i) = F*TComplex::Exp(phase); - xF.MakeReAndIm(i); + TComplex F(re_.at(i), im_.at(i)); + TComplex phase = TComplex(0., -TMath::TwoPi()*i*mu/(fSamples)); + xF.fCmplx.at(i) = F*TComplex::Exp(phase); + xF.MakeReAndIm(i); } - return; + return; } @@ -630,7 +800,7 @@ namespace opdet { * * Produce a waveform representing a guess of the input signal based on the * estimated nr of p.e. in the waveform peak. Based on the `fInputShape` - * variable, the input shape can be assumed as a δ-function scaled for + * variable, the input shape can be assumed as a δ-function scaled for * the nr of p.e. or as the LAr scintillation time profile * * @param s input signal @@ -638,12 +808,12 @@ namespace opdet { */ void Deconvolution::ComputeExpectedInput(std::vector& s, double nmax) { if (fInputShape == kScint) { - //ST profile + //ST profile auto const *Larprop = lar::providerFrom(); std::vectorSignalTime={(Larprop->ScintFastTimeConst()* 0.001),(Larprop->ScintSlowTimeConst()* 0.001)}; std::vectorSignalScint={Larprop->ScintYieldRatio(),1.-Larprop->ScintYieldRatio()}; - double dt = 1/fSampleFreq; + double dt = 1/fSampleFreq; double t = 0; for (size_t i=0; i()); // add a new empty waform + auto &spewfrm = fSinglePEWaveforms.back(); // get the reference to the waveform vector + std::string datafile; + // taking the file name as the first argument, + // the second argument is the local variable where to store the full path - both are std::string objects + sp.find_file(fname, datafile); + std::ifstream SPEData; + SPEData.open(datafile); + size_t n_columns = CountFileColumns(datafile.c_str()); + std::cout << "ncols= " << n_columns << std::endl; + if (fSPETemplateFileDataColumn >= n_columns) { + printf("Deconvolution::SourceSPETemplateFiles ERROR: "); + printf("The module is supposed to select column %lu, but only %lu columns are present.\n", + fSPETemplateFileDataColumn, n_columns); + throw art::Exception(art::errors::InvalidNumber); + } - Double_t buff[100] = {0}; + Double_t buff[100] = {0}; - std::string temp_str; - if (SPEData.is_open()) { - while (std::getline(SPEData, temp_str)) { - std::stringstream ss; ss << temp_str; - int icol = 0; - while (ss) {ss >> buff[icol]; ++icol;} + std::string temp_str; + if (SPEData.is_open()) { + while (std::getline(SPEData, temp_str)) { + std::stringstream ss; ss << temp_str; + int icol = 0; + while (ss) {ss >> buff[icol]; ++icol;} - fSinglePEWaveform.push_back(buff[fDigiDataColumn]); - } - } else { - printf("Deconvolution::produce ERROR "); - printf("Cannot open SPE template file.\n"); + spewfrm.push_back(buff[fSPETemplateFileDataColumn]); + } + } else { + printf("Deconvolution::produce ERROR "); + printf("Cannot open SPE template file.\n"); + + throw art::Exception(art::errors::FileOpenError); + } - throw art::Exception(art::errors::FileOpenError); + spewfrm.resize(fSamples, 0); + + SPEData.close(); + + // Set single p.e. maximum value + fSinglePEAmplitudes.push_back( TMath::Max(1.0, + *(std::max_element(spewfrm.begin(), spewfrm.end()))) ); + std::cout << "SPE Amplitude for template " << fSinglePEWaveforms.size() << ": " << fSinglePEAmplitudes.back() << std::endl; } + return; + } + + /** + * @brief Source template noise from files + * + * Source the noise templates from the dat files set by + * `fNoiseTemplateFiles`. + */ + void Deconvolution::SourceNoiseTemplateFiles() { + cet::search_path sp("FW_SEARCH_PATH"); + for (auto fname: fNoiseTemplateFiles) { + fNoiseTemplates.push_back(std::vector()); // add a new empty waform + auto &noisewfrm = fNoiseTemplates.back(); // get the reference to the waveform vector + std::string datafile; + // taking the file name as the first argument, + // the second argument is the local variable where to store the full path - both are std::string objects + sp.find_file(fname, datafile); + std::ifstream noiseData; + noiseData.open(datafile); + + std::string temp_str; + double temp = 0.; + if (noiseData.is_open()) { + while (std::getline(noiseData, temp_str)) { + std::stringstream ss; ss << temp_str; + ss >> temp; + noisewfrm.push_back(temp); + } + } else { + printf("Deconvolution::SourceNoiseTemplateFiles ERROR "); + printf("Cannot open noise template file.\n"); + + throw art::Exception(art::errors::FileOpenError); + } - fSinglePEWaveform.resize(fSamples, 0); + noisewfrm.resize(fSamples, 0.); - SPEData.close(); + noiseData.close(); - // Set single p.e. maximum value - fSinglePEAmplitude = TMath::Max(1.0, - *(std::max_element(fSinglePEWaveform.begin(), fSinglePEWaveform.end()))); - std::cout << "SPE Amplitude: " << fSinglePEAmplitude << std::endl; + } return; } + /** * @brief Count the nr of column in a txt file * @@ -730,54 +946,54 @@ namespace opdet { * @return nr of columns */ int Deconvolution::CountFileColumns(const char* file_path) { - std::ifstream file_; - file_.open(file_path); + std::ifstream file_; + file_.open(file_path); if (!file_.is_open()) { - printf("Deconvolution::CountFileColumns(%s) ERROR:\n", - file_path); - printf("Unable to open file."); - throw art::Exception(art::errors::FileOpenError); + printf("Deconvolution::CountFileColumns(%s) ERROR:\n", + file_path); + printf("Unable to open file."); + throw art::Exception(art::errors::FileOpenError); } - int N_COLUMNS = 0; - std::string line; - int iline = 0; + int N_COLUMNS = 0; + std::string line; + int iline = 0; while ( std::getline(file_, line) ) { - std::stringstream sstream; + std::stringstream sstream; sstream << line; std::string sub; - int n_columns = 0; + int n_columns = 0; while (sstream) { - sstream >> sub; - if (sub.length()) ++n_columns; + sstream >> sub; + if (sub.length()) ++n_columns; } if (iline == 1) {N_COLUMNS = n_columns;} else if (iline > 1) { if (n_columns != N_COLUMNS) { - printf("Deconvolution::CountFileColumns(%s): WARNING ", - file_path); - printf("Nr of columns change along the file!\n"); - N_COLUMNS = n_columns; + printf("Deconvolution::CountFileColumns(%s): WARNING ", + file_path); + printf("Nr of columns change along the file!\n"); + N_COLUMNS = n_columns; } } - iline++; + iline++; } - file_.close(); - return N_COLUMNS; + file_.close(); + return N_COLUMNS; } /** * @brief Compute normalization factor for a given filter * - * The filter normalization factor is obtained by applying the + * The filter normalization factor is obtained by applying the * filter to the single p.e. response template (i.e., to a 1 p.e. noiseless signal). - * the product of the filter should be a the best approximation of the - * input function achievable with the signal to noise ratio in such waveform, - * thus, the normalization consists in the integral of the product - * around the peak in a region defined as the one where the signal is positive. + * the product of the filter should be a the best approximation of the + * input function achievable with the signal to noise ratio in such waveform, + * thus, the normalization consists in the integral of the product + * around the peak in a region defined as the one where the signal is positive. * This factor is supposed to tend to 1 for high SNR (>= 10). * * @param xGH @@ -791,45 +1007,49 @@ namespace opdet { // time window const int shift = 0.5*fSamples; for (int i=0; iSetPointsComplex(&xGH.fRe[0], &xGH.fIm[0]); - fft_c2r->Transform(); - Double_t* x_ = fft_c2r->GetPointsReal(); + fft_c2r->SetPointsComplex(&xGH.fRe[0], &xGH.fIm[0]); + fft_c2r->Transform(); + Double_t* x_ = fft_c2r->GetPointsReal(); std::vector x(x_, x_ + fSamples); - Int_t imax = std::distance(x.begin(), std::max_element(x.begin(), x.end())); - Int_t ileft = imax; - Int_t iright = imax; + Int_t imax = std::distance(x.begin(), std::max_element(x.begin(), x.end())); + Int_t ileft = imax; + Int_t iright = imax; - while (x[ileft] > thrs && ileft > 0) ileft--; - while (x[iright] > thrs && iright < fSamples) iright++; + while (x[ileft] > thrs && ileft > 0) ileft--; + while (x[iright] > thrs && iright < fSamples) iright++; - double norm = 0; - for (Int_t k=ileft; k<=iright; k++) norm += x[k]; - norm /= (Double_t)fSamples; + double norm = 0; + for (Int_t k=ileft; k<=iright; k++) norm += x[k]; + norm /= (Double_t)fSamples; - if (norm > 1.0) norm = 1.0; + if (norm > 1.0) norm = 1.0; else if (norm <= 0.){ - printf("Deconvolution::ComputeAutoNormalization() WARNING: "); - printf(" bad normalization (%g), force to 1.0\n", norm); - norm = 1.0; + printf("Deconvolution::ComputeAutoNormalization() WARNING: "); + printf(" bad normalization (%g), force to 1.0\n", norm); + norm = 1.0; } - - return 1.0 / norm; + + return 1.0 / norm; } void Deconvolution::CopyToOutput(const std::vector& v, std::vector& target) { target.assign(v.begin(), v.end()); return; } + void Deconvolution::CopyToOutput(const std::vector& v, std::vector& target) { + target = std::vector(v.begin(), v.end()); + return; + } void Deconvolution::CopyToOutput(const CmplxWaveform_t& v, std::vector& target) { for (size_t i=0; i