diff --git a/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl b/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl index fecad004c..d9b166996 100755 --- a/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl +++ b/sbndcode/JobConfigurations/standard/reco/config/workflow_reco1.fcl @@ -18,18 +18,18 @@ sbnd_reco1_producers:{ rns: { module_type: RandomNumberSaver } ### optical deconvolution - opdecopmt: @local::SBNDOpDeconvolutionPMT - opdecoxarapuca: @local::SBNDOpDeconvolutionXARAPUCA + opdecopmt: @local::SBNDOpDeconvolutionPMT_realisticMC + opdecoxarapuca: @local::SBNDOpDeconvolutionXARAPUCA ### optical hit finders # ophit: @local::sbnd_hit_finder - ophitpmt: @local::SBNDDecoOpHitFinderPMT - ophitxarapuca: @local::SBNDDecoOpHitFinderXArapuca + ophitpmt: @local::SBNDDecoOpHitFinderPMT_realisticMC + ophitxarapuca: @local::SBNDDecoOpHitFinderXArapuca ### flash finders # opflash: @local::sbnd_opflash - opflashtpc0: @local::SBNDDecoSimpleFlashTPC0 - opflashtpc1: @local::SBNDDecoSimpleFlashTPC1 + opflashtpc0: @local::SBNDDecoSimpleFlashTPC0_realisticMC + opflashtpc1: @local::SBNDDecoSimpleFlashTPC1_realisticMC ## opflash(arapucas): @local::sbnd_opflash opflashtpc0xarapuca: @local::SBNDDecoSimpleFlashTPC0Arapuca diff --git a/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl b/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl index e5552ebc4..4b6ae0e4d 100755 --- a/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl +++ b/sbndcode/JobConfigurations/standard/reco/reco1_data.fcl @@ -1,8 +1,8 @@ #include "drops_reco1.fcl" #include "wcsp_data_sbnd.fcl" -#include "opdeconvolution_sbnd_data.fcl" -#include "sbnd_ophitfinder_deco_data.fcl" -#include "sbnd_flashfinder_deco_data.fcl" +#include "opdeconvolution_sbnd.fcl" +#include "sbnd_ophitfinder_deco.fcl" +#include "sbnd_flashfinder_deco.fcl" #include "wfalign_sbnd_data.fcl" #include "standard_reco1_sbnd.fcl" diff --git a/sbndcode/JobConfigurations/standard/standard_detsim_sbnd.fcl b/sbndcode/JobConfigurations/standard/standard_detsim_sbnd.fcl index cd9094870..2dabe25f4 100755 --- a/sbndcode/JobConfigurations/standard/standard_detsim_sbnd.fcl +++ b/sbndcode/JobConfigurations/standard/standard_detsim_sbnd.fcl @@ -34,6 +34,7 @@ #include "wcsimsp_sbnd.fcl" #include "crtsimmodules_sbnd.fcl" #include "opdetdigitizer_sbnd.fcl" +#include "sbnd_pmtpulseoscillation_config.fcl" #include "rootoutput_sbnd.fcl" #include "pmtmcmetricproducer.fcl" @@ -63,15 +64,16 @@ physics: producers: { - rns: { module_type: "RandomNumberSaver" } - simtpc2d: @local::sbnd_wcls_simsp_bothrois - crtsim: @local::sbnd_crtsim - opdaq: @local::sbnd_opdetdigitizer + rns: { module_type: "RandomNumberSaver" } + simtpc2d: @local::sbnd_wcls_simsp_bothrois + crtsim: @local::sbnd_crtsim + opdaq: @local::sbnd_opdetdigitizer pmtmetricmc: @local::sbnd_pmtmcmetricproducer + pmtpulseoscillation: @local::sbnd_pmtpulseoscillation } #define the producer and filter modules for this path, order matters, - simulate: [rns, simtpc2d, crtsim, opdaq, pmtmetricmc] + simulate: [rns, simtpc2d, crtsim, opdaq, pmtmetricmc, pmtpulseoscillation] #define the output stream, there could be more than one if using filters stream1: [ out1 ] diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd.fcl index 9599a707c..0d7790bf5 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd.fcl +++ b/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd.fcl @@ -1,4 +1,5 @@ #include "opdeconvolution_alg.fcl" +#include "opdeconvolution_alg_data.fcl" BEGIN_PROLOG @@ -11,6 +12,8 @@ SBNDOpDeconvolution: OpDecoAlg: @local::OpDeconvolutionAlg } +###### PMT IDEAL MC ###### + SBNDOpDeconvolutionPMT: @local::SBNDOpDeconvolution SBNDOpDeconvolutionPMT.PDTypes: ["pmt_coated", "pmt_uncoated"] SBNDOpDeconvolutionPMT.Electronics: [""] @@ -20,6 +23,34 @@ SBNDOpDeconvolutionPMT.OpDecoAlg.FilterParams: [0.049, 2] #Freq in GHz SBNDOpDeconvolutionPMT.OpDecoAlg.Filter: "(x>0)*exp(-0.5*pow(x/[0],[1]))" #Gauss filter, remove DC component F(0)=0 SBNDOpDeconvolutionPMT.OpDecoAlg.DecoWaveformPrecision: 0.005 + +###### PMT DATA ###### + +SBNDOpDeconvolutionPMT_data: @local::SBNDOpDeconvolution +SBNDOpDeconvolutionPMT_data.OpDecoAlg: @local::OpDeconvolutionAlgData +SBNDOpDeconvolutionPMT_data.InputLabel: "wfalign:PMTChannels" +SBNDOpDeconvolutionPMT_data.PDTypes: ["pmt_coated", "pmt_uncoated"] +SBNDOpDeconvolutionPMT_data.Electronics: [""] +SBNDOpDeconvolutionPMT_data.OpDecoAlg.OpDetDataFile: "./OpDetSim/digi_pmt_sbnd_data_OV6.root" +SBNDOpDeconvolutionPMT_data.OpDecoAlg.Filter: "(x>0)*exp(-0.5*pow(x/[0],[1]))" #Gauss filter, remove DC component F(0)=0 +SBNDOpDeconvolutionPMT_data.OpDecoAlg.DecoWaveformPrecision: 0.005 +SBNDOpDeconvolutionPMT_data.OpDecoAlg.SkipChannelList: [39, 66, 67, 71, 85, 86, 87, 92, 115, 138, 141, 170, 197, 217, 218, 221, 222, 223, 226, 245, 248, 249, 302] + +###### PMT REALISTIC MC ###### + +SBNDOpDeconvolutionPMT_realisticMC: @local::SBNDOpDeconvolution +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg: @local::OpDeconvolutionAlgData +SBNDOpDeconvolutionPMT_realisticMC.InputLabel: "pmtpulseoscillation" +SBNDOpDeconvolutionPMT_realisticMC.PDTypes: ["pmt_coated", "pmt_uncoated"] +SBNDOpDeconvolutionPMT_realisticMC.Electronics: [""] +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg.OpDetDataFile: "./OpDetSim/digi_pmt_sbnd_data_OV6.root" +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg.Filter: "(x>0)*exp(-0.5*pow(x/[0],[1]))" #Gauss filter, remove DC component F(0)=0 +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg.DecoWaveformPrecision: 0.005 +SBNDOpDeconvolutionPMT_realisticMC.OpDecoAlg.SkipChannelList: [39, 66, 67, 71, 85, 86, 87, 92, 115, 138, 141, 170, 197, 217, 218, 221, 222, 223, 226, 245, 248, 249, 302] + + +###### XA IDEAL MC ###### + SBNDOpDeconvolutionXARAPUCA: @local::SBNDOpDeconvolution SBNDOpDeconvolutionXARAPUCA.PDTypes: ["xarapuca_vuv", "xarapuca_vis"] SBNDOpDeconvolutionXARAPUCA.Electronics: ["daphne"] diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd_data.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd_data.fcl deleted file mode 100644 index 66ad23bc2..000000000 --- a/sbndcode/OpDetReco/OpDeconvolution/job/opdeconvolution_sbnd_data.fcl +++ /dev/null @@ -1,24 +0,0 @@ -#include "opdeconvolution_alg_data.fcl" - -BEGIN_PROLOG - -SBNDOpDeconvolution: -{ - module_type: "SBNDOpDeconvolution" - InputLabel: "wfalign:PMTChannels" - PDTypes: [] - Electronics: [] - OpDecoAlg: @local::OpDeconvolutionAlgData -} - -SBNDOpDeconvolutionPMT_data: @local::SBNDOpDeconvolution -SBNDOpDeconvolutionPMT_data.PDTypes: ["pmt_coated", "pmt_uncoated"] -SBNDOpDeconvolutionPMT_data.Electronics: [""] -SBNDOpDeconvolutionPMT_data.OpDecoAlg.OpDetDataFile: "./OpDetSim/digi_pmt_sbnd_data_OV6.root" -#SBNDOpDeconvolutionPMT_data.OpDecoAlg.UseParamFilter: true -#SBNDOpDeconvolutionPMT_data.OpDecoAlg.FilterParams: [0.049, 2] #Freq in GHz -SBNDOpDeconvolutionPMT_data.OpDecoAlg.Filter: "(x>0)*exp(-0.5*pow(x/[0],[1]))" #Gauss filter, remove DC component F(0)=0 -SBNDOpDeconvolutionPMT_data.OpDecoAlg.DecoWaveformPrecision: 0.005 -SBNDOpDeconvolutionPMT_data.OpDecoAlg.SkipChannelList: [39, 66, 67, 71, 85, 86, 87, 92, 115, 138, 141, 170, 197, 217, 218, 221, 222, 223, 226, 245, 248, 249, 302] - -END_PROLOG diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder_data.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder_data.fcl index 9d8c2af9d..a2ebd0b4a 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder_data.fcl +++ b/sbndcode/OpDetReco/OpDeconvolution/job/run_decohitfinder_data.fcl @@ -3,11 +3,10 @@ #include "messages_sbnd.fcl" #include "sam_sbnd.fcl" -#include "opdeconvolution_sbnd_data.fcl" -#include "sbnd_flashfinder_deco_data.fcl" -#include "sbnd_ophitfinder_deco_data.fcl" +#include "opdeconvolution_sbnd.fcl" +#include "sbnd_flashfinder_deco.fcl" +#include "sbnd_ophitfinder_deco.fcl" #include "wfalign_sbnd_data.fcl" - process_name: DecoRecoData services: @@ -57,7 +56,6 @@ physics: opflashtpc1 ] - stream1: [out1] trigger_paths: [reco] diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco.fcl index fd231e3f8..41d453e49 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco.fcl +++ b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco.fcl @@ -2,7 +2,11 @@ BEGIN_PROLOG + ####OpFlash finder for PMT deconvolved waveforms##### + + +###### PMT IDEAL MC ###### ###TPC0 SBNDDecoSimpleFlashTPC0: @local::SBNDSimpleFlashTPC0 SBNDDecoSimpleFlashTPC0.PECalib.SPEAreaGain: 200 @@ -22,8 +26,57 @@ SBNDDecoSimpleFlashTPC1.ReadoutDelay: 0.135 //cable time delay in us SBNDDecoSimpleFlashTPC1.CorrectLightPropagation: true +###### PMT DATA ###### +###TPC0 +SBNDDecoSimpleFlashTPC0_data: @local::SBNDSimpleFlashTPC0 +SBNDDecoSimpleFlashTPC0_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" +SBNDDecoSimpleFlashTPC0_data.DriftEstimatorConfig.DataCalibration: true +SBNDDecoSimpleFlashTPC0_data.PECalib.SPEAreaGain: 200 +SBNDDecoSimpleFlashTPC0_data.OpHitProducers: ["ophitpmt"] +SBNDDecoSimpleFlashTPC0_data.OpHitInputTime: "RiseTime" +SBNDDecoSimpleFlashTPC0_data.UseT0Tool: true +SBNDDecoSimpleFlashTPC0_data.ReadoutDelay: 0 //cable time delay in us +SBNDDecoSimpleFlashTPC0_data.CorrectLightPropagation: true +SBNDDecoSimpleFlashTPC0_data.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data_v2.root" + +#TPC1 +SBNDDecoSimpleFlashTPC1_data: @local::SBNDSimpleFlashTPC1 +SBNDDecoSimpleFlashTPC1_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" +SBNDDecoSimpleFlashTPC1_data.DriftEstimatorConfig.DataCalibration: true +SBNDDecoSimpleFlashTPC1_data.PECalib.SPEAreaGain: 200 +SBNDDecoSimpleFlashTPC1_data.OpHitProducers: ["ophitpmt"] +SBNDDecoSimpleFlashTPC1_data.OpHitInputTime: "RiseTime" +SBNDDecoSimpleFlashTPC1_data.UseT0Tool: true +SBNDDecoSimpleFlashTPC1_data.ReadoutDelay: 0 //cable time delay in us +SBNDDecoSimpleFlashTPC1_data.CorrectLightPropagation: true + + +###### PMT REALISTIC MC ###### +###TPC0 +SBNDDecoSimpleFlashTPC0_realisticMC: @local::SBNDSimpleFlashTPC0 +SBNDDecoSimpleFlashTPC0_realisticMC.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" +SBNDDecoSimpleFlashTPC0_realisticMC.DriftEstimatorConfig.DataCalibration: true +SBNDDecoSimpleFlashTPC0_realisticMC.PECalib.SPEAreaGain: 200 +SBNDDecoSimpleFlashTPC0_realisticMC.OpHitProducers: ["ophitpmt"] +SBNDDecoSimpleFlashTPC0_realisticMC.OpHitInputTime: "RiseTime" +SBNDDecoSimpleFlashTPC0_realisticMC.UseT0Tool: true +SBNDDecoSimpleFlashTPC0_realisticMC.ReadoutDelay: 0 //cable time delay in us +SBNDDecoSimpleFlashTPC0_realisticMC.CorrectLightPropagation: true + +#TPC1 +SBNDDecoSimpleFlashTPC1_realisticMC: @local::SBNDSimpleFlashTPC1 +SBNDDecoSimpleFlashTPC1_realisticMC.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" +SBNDDecoSimpleFlashTPC1_realisticMC.DriftEstimatorConfig.DataCalibration: true +SBNDDecoSimpleFlashTPC1_realisticMC.PECalib.SPEAreaGain: 200 +SBNDDecoSimpleFlashTPC1_realisticMC.OpHitProducers: ["ophitpmt"] +SBNDDecoSimpleFlashTPC1_realisticMC.OpHitInputTime: "RiseTime" +SBNDDecoSimpleFlashTPC1_realisticMC.UseT0Tool: true +SBNDDecoSimpleFlashTPC1_realisticMC.ReadoutDelay: 0 //cable time delay in us +SBNDDecoSimpleFlashTPC1_realisticMC.CorrectLightPropagation: true +SBNDDecoSimpleFlashTPC1_realisticMC.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data_v2.root" +###### XA IDEAL MC ###### ####OpFlash finder for XArapucas deconvolved waveforms##### ###TPC0 diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco_data.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco_data.fcl deleted file mode 100644 index 45b3aa338..000000000 --- a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_flashfinder_deco_data.fcl +++ /dev/null @@ -1,27 +0,0 @@ -#include "sbnd_flashfinder_data.fcl" - -BEGIN_PROLOG - -####OpFlash finder for PMT deconvolved waveforms##### -###TPC0 -SBNDDecoSimpleFlashTPC0_data: @local::SBNDSimpleFlashTPC0_data -SBNDDecoSimpleFlashTPC0_data.PECalib.SPEAreaGain: 200 -SBNDDecoSimpleFlashTPC0_data.OpHitProducers: ["ophitpmt"] -SBNDDecoSimpleFlashTPC0_data.OpHitInputTime: "RiseTime" -SBNDDecoSimpleFlashTPC0_data.UseT0Tool: true -SBNDDecoSimpleFlashTPC0_data.ReadoutDelay: 0 //cable time delay in us -SBNDDecoSimpleFlashTPC0_data.CorrectLightPropagation: true -SBNDDecoSimpleFlashTPC0_data.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data.root" - -#TPC1 -SBNDDecoSimpleFlashTPC1_data: @local::SBNDSimpleFlashTPC1_data -SBNDDecoSimpleFlashTPC1_data.PECalib.SPEAreaGain: 200 -SBNDDecoSimpleFlashTPC1_data.OpHitProducers: ["ophitpmt"] -SBNDDecoSimpleFlashTPC1_data.OpHitInputTime: "RiseTime" -SBNDDecoSimpleFlashTPC1_data.UseT0Tool: true -SBNDDecoSimpleFlashTPC1_data.ReadoutDelay: 0 //cable time delay in us -SBNDDecoSimpleFlashTPC1_data.CorrectLightPropagation: true -SBNDDecoSimpleFlashTPC1_data.DriftEstimatorConfig.CalibrationFile: "OpDetReco/PMTRatioCalibration_data.root" - - -END_PROLOG diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco.fcl index 2911860a9..989fe0193 100644 --- a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco.fcl +++ b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco.fcl @@ -2,6 +2,9 @@ BEGIN_PROLOG +###### MC CONFIGURATION ###### + + #####OpHit finder for PMT deconvolved waveforms##### SBNDDecoOpHitFinderPMT: @local::sbnd_ophit_finder_pmt ####SPE area must be 1./DecoWaveformPrecision @@ -63,4 +66,44 @@ SBNDDecoOpHitFinderXArapuca.PedAlgoPset.NumSampleTail: 50 SBNDDecoOpHitFinderXArapuca.PedAlgoPset.Method:2 +###### DATA CONFIGURATION ###### + +#####OpHit finder for PMT deconvolved waveforms##### +SBNDDecoOpHitFinderPMT_data: @local::sbnd_ophit_finder_pmt +####SPE area must be 1./DecoWaveformPrecision +SBNDDecoOpHitFinderPMT_data.SPEArea: 200 +SBNDDecoOpHitFinderPMT_data.InputModule: "opdecopmt" +SBNDDecoOpHitFinderPMT_data.HitThreshold: 1 +SBNDDecoOpHitFinderPMT_data.RiseTimeCalculator: @local::sbnd_opreco_risetimecalculator + +#HitAlgoPset +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThreshold: 25 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThresholdByChannel: true +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThresholdVector: [50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 18.0, 35.0, 19.0, 21.0, 19.0, 20.0, 21.0, 27.0, 19.0, 27.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 19.0, 37.0, 50.0, 20.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 24.0, 19.0, 26.0, 19.0, 12.0, 30.0, 50.0, 25.0, 14.0, 12.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 50.0, 50.0, 50.0, 19.0, 19.0, 24.0, 21.0, 28.0, 19.0, 22.0, 18.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 50.0, 19.0, 22.0, 35.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 30.0, 50.0, 20.0, 25.0, 23.0, 32.0, 20.0, 22.0, 20.0, 18.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 23.0, 50.0, 19.0, 25.0, 19.0, 18.0, 19.0, 15.0, 21.0, 24.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 20.0, 30.0, 26.0, 18.0, 34.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 27.0, 50.0, 50.0, 21.0, 28.0, 50.0, 50.0, 50.0, 19.0, 19.0, 50.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 19.0, 30.0, 20.0, 18.0, 50.0, 18.0, 19.0, 13.0, 50.0, 25.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 25.0, 22.0, 19.0, 29.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 23.0, 29.0, 19.0, 27.0, 20.0, 25.0, 22.0, 50.0, 19.0, 25.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0] +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NSigmaThreshold: 3.4 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.EndADCThreshold: 8 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.EndNSigmaThreshold: 0.2 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.MinPulseWidth: 4 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.Name: "SlidingWindow" +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NumPostSample: 6 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NumPreSample: 3 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.PositivePolarity: true +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.TailADCThreshold: 2 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.TailNSigmaThreshold: 2 +SBNDDecoOpHitFinderPMT_data.HitAlgoPset.Verbosity: false + +#BaselinePset +SBNDDecoOpHitFinderPMT_data.PedAlgoPset.Name:"Edges" +SBNDDecoOpHitFinderPMT_data.PedAlgoPset.NumSampleFront:200 +SBNDDecoOpHitFinderPMT_data.PedAlgoPset.NumSampleTail:200 +SBNDDecoOpHitFinderPMT_data.PedAlgoPset.Method:2 + + +###### REALISTICMC CONFIGURATION ###### + +#####OpHit finder for PMT deconvolved waveforms##### +SBNDDecoOpHitFinderPMT_realisticMC: @local::SBNDDecoOpHitFinderPMT_data +####SPE area must be 1./DecoWaveformPrecision +SBNDDecoOpHitFinderPMT_realisticMC.InputModule: "pmtpulseoscillation" + END_PROLOG diff --git a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco_data.fcl b/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco_data.fcl deleted file mode 100644 index 02eadf426..000000000 --- a/sbndcode/OpDetReco/OpDeconvolution/job/sbnd_ophitfinder_deco_data.fcl +++ /dev/null @@ -1,36 +0,0 @@ -#include "ophit_finder_sbnd.fcl" - -BEGIN_PROLOG - -#####OpHit finder for PMT deconvolved waveforms##### -SBNDDecoOpHitFinderPMT_data: @local::sbnd_ophit_finder_pmt -####SPE area must be 1./DecoWaveformPrecision -SBNDDecoOpHitFinderPMT_data.SPEArea: 200 -SBNDDecoOpHitFinderPMT_data.InputModule: "opdecopmt" -SBNDDecoOpHitFinderPMT_data.HitThreshold: 1 -SBNDDecoOpHitFinderPMT_data.RiseTimeCalculator: @local::sbnd_opreco_risetimecalculator - -#HitAlgoPset -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThreshold: 25 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThresholdByChannel: true -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.ADCThresholdVector: [50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 18.0, 35.0, 19.0, 21.0, 19.0, 20.0, 21.0, 27.0, 19.0, 27.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 19.0, 37.0, 50.0, 20.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 24.0, 19.0, 26.0, 19.0, 12.0, 30.0, 50.0, 25.0, 14.0, 12.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 50.0, 50.0, 50.0, 19.0, 19.0, 24.0, 21.0, 28.0, 19.0, 22.0, 18.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 50.0, 19.0, 22.0, 35.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 30.0, 50.0, 20.0, 25.0, 23.0, 32.0, 20.0, 22.0, 20.0, 18.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 23.0, 50.0, 19.0, 25.0, 19.0, 18.0, 19.0, 15.0, 21.0, 24.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 20.0, 30.0, 26.0, 18.0, 34.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 27.0, 50.0, 50.0, 21.0, 28.0, 50.0, 50.0, 50.0, 19.0, 19.0, 50.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 19.0, 30.0, 20.0, 18.0, 50.0, 18.0, 19.0, 13.0, 50.0, 25.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 19.0, 25.0, 22.0, 19.0, 29.0, 19.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 22.0, 23.0, 29.0, 19.0, 27.0, 20.0, 25.0, 22.0, 50.0, 19.0, 25.0, 21.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0] -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NSigmaThreshold: 3.4 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.EndADCThreshold: 8 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.EndNSigmaThreshold: 0.2 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.MinPulseWidth: 4 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.Name: "SlidingWindow" -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NumPostSample: 6 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.NumPreSample: 3 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.PositivePolarity: true -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.TailADCThreshold: 2 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.TailNSigmaThreshold: 2 -SBNDDecoOpHitFinderPMT_data.HitAlgoPset.Verbosity: false - -#BaselinePset -SBNDDecoOpHitFinderPMT_data.PedAlgoPset.Name:"Edges" -SBNDDecoOpHitFinderPMT_data.PedAlgoPset.NumSampleFront:200 -SBNDDecoOpHitFinderPMT_data.PedAlgoPset.NumSampleTail:200 -SBNDDecoOpHitFinderPMT_data.PedAlgoPset.Method:2 - - -END_PROLOG diff --git a/sbndcode/OpDetReco/OpFlash/job/sbnd_flashfinder_data.fcl b/sbndcode/OpDetReco/OpFlash/job/sbnd_flashfinder_data.fcl deleted file mode 100644 index 19c4ccfa1..000000000 --- a/sbndcode/OpDetReco/OpFlash/job/sbnd_flashfinder_data.fcl +++ /dev/null @@ -1,36 +0,0 @@ -#include "sbnd_flashalgo.fcl" -#include "sbnd_flashcalib.fcl" -#include "sbnd_flashgeoalgo.fcl" -#include "sbnd_flasht0algo.fcl" -#include "sbnd_driftestimatoralgo.fcl" - -BEGIN_PROLOG - -SBNDSimpleFlash: -{ - module_type : "SBNDFlashFinder" - FlashFinderAlgo : "SimpleFlashAlgo" - AlgoConfig : @local::SimpleFlashStandard - OpHitProducers : ["ophitpmt","ophitarapuca"] - OpFlashProducer : "opflash" - PECalib : @local::NoCalib - OpHitInputTime : "PeakTime" - FlashGeoConfig : @local::FlashGeoThreshold - UseT0Tool : false - FlashT0Config : @local::FlashT0SelectedChannels - ReadoutDelay : 0. //in us - CorrectLightPropagation : false - DriftEstimatorConfig : @local::DriftEstimatorPMTRatio -} - -SBNDSimpleFlashTPC0_data: @local::SBNDSimpleFlash -SBNDSimpleFlashTPC0_data.AlgoConfig: @local::SimpleFlashTPC0 -SBNDSimpleFlashTPC0_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" -SBNDSimpleFlashTPC0_data.DriftEstimatorConfig.DataCalibration: true - -SBNDSimpleFlashTPC1_data: @local::SBNDSimpleFlash -SBNDSimpleFlashTPC1_data.AlgoConfig: @local::SimpleFlashTPC1 -SBNDSimpleFlashTPC1_data.DriftEstimatorConfig.tool_type: "DriftEstimatorPMTRatio" -SBNDSimpleFlashTPC1_data.DriftEstimatorConfig.DataCalibration: true - -END_PROLOG diff --git a/sbndcode/OpDetSim/CMakeLists.txt b/sbndcode/OpDetSim/CMakeLists.txt index 796328b68..35c22cf77 100755 --- a/sbndcode/OpDetSim/CMakeLists.txt +++ b/sbndcode/OpDetSim/CMakeLists.txt @@ -177,6 +177,7 @@ install_source() cet_enable_asserts() add_subdirectory(PMTAlg) add_subdirectory(HDWvf) +add_subdirectory(PMTPulseOscillations) # install sbnd_pds_mapping.json with mapping of the photon detectors install_fw(LIST sbnd_pds_mapping.json) diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc index 53ec18074..345a80bb2 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.cc @@ -17,6 +17,8 @@ namespace opdet { , fPMTCoatedVISEff_tpc1(fParams.PMTCoatedVISEff_tpc1 / fParams.larProp->ScintPreScale()) , fPMTUncoatedEff_tpc1(fParams.PMTUncoatedEff_tpc1/ fParams.larProp->ScintPreScale()) // , fSinglePEmodel(fParams.SinglePEmodel) + , fUseDataNoise(fParams.UseDataNoise) + , fOpDetNoiseFile(fParams.OpDetNoiseFile) , fEngine(fParams.engine) , fFlatGen(*fEngine) , fPoissonQGen(*fEngine) @@ -45,7 +47,6 @@ namespace opdet { fSampling = fSampling / 1000.0; //in GHz, to cancel with ns fSamplingPeriod = 1./fSampling; - std::string fname; cet::search_path sp("FW_SEARCH_PATH"); sp.find_file(fParams.PMTDataFile, fname); @@ -56,54 +57,117 @@ namespace opdet { file->GetObject("timeTPB", timeTPB_p); fTimeTPB = std::make_unique (*fEngine, timeTPB_p->data(), timeTPB_p->size()); + + // PMT calibration database service + fPMTCalibrationDatabaseService = lar::providerFrom(); + fPMTHDOpticalWaveformsPtr = art::make_tool(fParams.HDOpticalWaveformParams); + int NOpChannels = fWireReadout.NOpChannels(); + //Resize the SER vector to the number of channels + fSinglePEWave_HD.resize(NOpChannels); //shape of single pulse - if (fParams.PMTSinglePEmodel) { + if (fParams.PMTSinglePEmodel=="testbench") { mf::LogDebug("DigiPMTSBNDAlg") << " using testbench pe response"; std::vector* SinglePEVec_p; file->GetObject("SinglePEVec_HD", SinglePEVec_p); fSinglePEWave = *SinglePEVec_p; - // Prepare HD waveforms - fPMTHDOpticalWaveformsPtr = art::make_tool(fParams.HDOpticalWaveformParams); - fPMTHDOpticalWaveformsPtr->produceSER_HD(fSinglePEWave_HD,fSinglePEWave); - - pulsesize = fSinglePEWave_HD[0].size(); + // Create the same SER for all channels + for(size_t i=0; iproduceSER_HD(fSinglePEWave_HD[i], fSinglePEWave); + } + pulsesize = fSinglePEWave_HD[0][0].size(); mf::LogDebug("DigiPMTSBNDAlg")<<"HD wvfs size: "<getReconstructChannel(i)){ + dataSERSize = fPMTCalibrationDatabaseService->getSER(i).size(); + break; + } + } + fAverageDataSER.resize(dataSERSize); + fAverageDataSPEAmplitude = 0; + int nChannelsWithDataSER = 0; + for(size_t i=0; igetReconstructChannel(i)) continue; + for (size_t j = 0; j < fAverageDataSER.size() && j < fPMTCalibrationDatabaseService->getSER(i).size(); ++j) { fAverageDataSER[j] += fPMTCalibrationDatabaseService->getSER(i)[j]; } + fAverageDataSPEAmplitude += fPMTCalibrationDatabaseService->getSPEAmplitude(i); + nChannelsWithDataSER++; + } + std::transform(fAverageDataSER.begin(), fAverageDataSER.end(), fAverageDataSER.begin(), [nChannelsWithDataSER](double val) {return val / nChannelsWithDataSER;}); + double PeakValue = *std::max_element(fAverageDataSER.begin(), fAverageDataSER.end(), [](double a, double b) {return std::abs(a) < std::abs(b);}); + fAverageDataSPEAmplitude = fAverageDataSPEAmplitude / nChannelsWithDataSER; + // Normalise the average data SER to the average SPE amplitude + double AverageDataPENormalization = std::abs(fAverageDataSPEAmplitude/PeakValue); + std::transform(fAverageDataSER.begin(), fAverageDataSER.end(), fAverageDataSER.begin(), [AverageDataPENormalization](double val) {return val * AverageDataPENormalization;}); + // Read the SER from the calibration database (channel independent) + for(size_t i=0; igetOnPMT(i)){ + // PMT is OFF + continue; + } + else if( fPMTCalibrationDatabaseService->getOnPMT(i) && !fPMTCalibrationDatabaseService->getReconstructChannel(i)){ + // PMT is ON but does not have a calibrated SER. These are only relevant for trigger emulation, not used in downstream reconstruciton + fPMTHDOpticalWaveformsPtr->produceSER_HD(fSinglePEWave_HD[i], fAverageDataSER); + } + else{ + // PMT is ON and has a calibrated SER + fSinglePEWave = fPMTCalibrationDatabaseService->getSER(i); + double SPEAmplitude = fPMTCalibrationDatabaseService->getSPEAmplitude(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(SPEAmplitude/SPEPeakValue); + std::transform(fSinglePEWave.begin(), fSinglePEWave.end(), fSinglePEWave.begin(), [SinglePENormalization](double val) {return val * SinglePENormalization;}); + if(fSinglePEWave.size()==0) continue; + if(fSinglePEWave.size()>0) pulsesize = fSinglePEWave.size(); + fPMTHDOpticalWaveformsPtr->produceSER_HD(fSinglePEWave_HD[i], fSinglePEWave); + } + } + } + else{ + throw cet::exception("DigiPMTSBNDAlg") << "Wrong PMTSinglePEmodel configuration: " << fParams.PMTSinglePEmodel << std::endl; + } if(fParams.MakeGainFluctuations){ fPMTGainFluctuationsPtr = art::make_tool(fParams.GainFluctuationsParams); } if(fParams.SimulateNonLinearity){ - fPMTNonLinearityPtr = art::make_tool(fParams.NonLinearityParams); - } - + fPMTNonLinearityPtr = art::make_tool(fParams.NonLinearityParams); + fPMTNonLinearityPtr->ConfigureNonLinearity(); + } // infer pulse polarity from SER peak sign double minADC_SinglePE = *min_element(fSinglePEWave.begin(), fSinglePEWave.end()); double maxADC_SinglePE = *max_element(fSinglePEWave.begin(), fSinglePEWave.end()); - fPositivePolarity = std::abs(maxADC_SinglePE) > std::abs(minADC_SinglePE); - + fPositivePolarity = std::abs(maxADC_SinglePE) > std::abs(minADC_SinglePE); + // get ADC saturation value // currently assumes all dynamic range for PE (no overshoot) fADCSaturation = (fPositivePolarity ? fParams.PMTBaseline + fParams.PMTADCDynamicRange : fParams.PMTBaseline - fParams.PMTADCDynamicRange); - - file->Close(); + // Initialize noise file + std::string fname_noise; + cet::search_path sp_noise("FW_SEARCH_PATH"); + if(fParams.UseDataNoise){ + sp_noise.find_file(fParams.OpDetNoiseFile, fname_noise); + noise_file = TFile::Open(fname_noise.c_str(), "READ"); + } } // end constructor - - DigiPMTSBNDAlg::~DigiPMTSBNDAlg(){} + DigiPMTSBNDAlg::~DigiPMTSBNDAlg(){ + if(fParams.UseDataNoise) noise_file->Close(); + } void DigiPMTSBNDAlg::ConstructWaveformUncoatedPMT( @@ -186,7 +250,7 @@ namespace opdet { //photon time in ns (w.r.t. the waveform start time a.k.a t_min) tphoton = ttsTime + simphotons[i].Time - t_min + ttpb + fParams.CableTime; - // store the pgoton time if it's within the readout window + // store the photon time if it's within the readout window if(tphoton > 0 && tphoton < nPE_v.size()) nPE_v[(size_t)tphoton]++; } } @@ -194,16 +258,22 @@ namespace opdet { for(size_t t=0; t 0) { if(fParams.SimulateNonLinearity){ - AddSPE(t, wave, fPMTNonLinearityPtr->NObservedPE(t, nPE_v) ); + AddSPE(t, wave, ch, fPMTNonLinearityPtr->NObservedPE(ch, t, nPE_v) ); } else{ - AddSPE(t, wave, nPE_v[t]); + AddSPE(t, wave, ch, nPE_v[t]); } } } - if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave); - if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave); + if(fParams.UseDataNoise) AddDataNoise(wave, ch); + else + { + if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave, ch); + } + + if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave, ch); + CreateSaturation(wave); } @@ -265,17 +335,22 @@ namespace opdet { for(size_t t=0; t 0) { if(fParams.SimulateNonLinearity){ - AddSPE(t, wave, fPMTNonLinearityPtr->NObservedPE(t, nPE_v) ); + AddSPE(t, wave, ch, fPMTNonLinearityPtr->NObservedPE(ch, t, nPE_v) ); } else{ - AddSPE(t, wave, nPE_v[t]); + AddSPE(t, wave, ch, nPE_v[t]); } } } //Adding noise and saturation - if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave); - if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave); + if(fParams.UseDataNoise) AddDataNoise(wave, ch); + else + { + if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave, ch); + } + if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave, ch); + CreateSaturation(wave); } @@ -320,16 +395,21 @@ namespace opdet { for(size_t t=0; t 0) { if(fParams.SimulateNonLinearity){ - AddSPE(t, wave, fPMTNonLinearityPtr->NObservedPE(t, nPE_v) ); + AddSPE(t, wave, ch, fPMTNonLinearityPtr->NObservedPE(ch, t, nPE_v) ); } else{ - AddSPE(t, wave, nPE_v[t]); + AddSPE(t, wave, ch, nPE_v[t]); } } } - if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave); - if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave); + if(fParams.UseDataNoise) AddDataNoise(wave, ch); + else + { + if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave, ch); + } + if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave, ch); + CreateSaturation(wave); } @@ -398,17 +478,22 @@ namespace opdet { for(size_t t=0; t 0) { if(fParams.SimulateNonLinearity){ - AddSPE(t, wave, fPMTNonLinearityPtr->NObservedPE(t, nPE_v) ); + AddSPE(t, wave, ch, fPMTNonLinearityPtr->NObservedPE(ch, t, nPE_v) ); } else{ - AddSPE(t, wave, nPE_v[t]); + AddSPE(t, wave, ch, nPE_v[t]); } } } //Adding noise and saturation - if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave); - if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave); + if(fParams.UseDataNoise) AddDataNoise(wave, ch); + else + { + if(fParams.PMTBaselineRMS > 0.0) AddLineNoise(wave, ch); + } + if(fParams.PMTDarkNoiseRate > 0.0) AddDarkNoise(wave, ch); + CreateSaturation(wave); } @@ -438,8 +523,9 @@ namespace opdet { } - void DigiPMTSBNDAlg::AddSPE(size_t time, std::vector& wave, double npe) - { + void DigiPMTSBNDAlg::AddSPE(size_t time, std::vector& wave, int ch, double npe) + { + if(!fPMTCalibrationDatabaseService->getOnPMT(ch)) return; // time bin HD (double precision) // used to gert the time-shifted SER double time_bin_hd = fSampling*time; @@ -454,11 +540,11 @@ namespace opdet { // simulate gain fluctuations double npe_anode = npe; if(fParams.MakeGainFluctuations) - npe_anode=fPMTGainFluctuationsPtr->GainFluctuation(npe, fEngine); - + //npe_anode=fPMTGainFluctuationsPtr->GainFluctuation(npe, fEngine); + npe_anode=fPMTGainFluctuationsPtr->GainFluctuation(ch, npe, fEngine); // add SER to the waveform std::transform(min_it, max_it, - fSinglePEWave_HD[wvf_shift].begin(), min_it, + fSinglePEWave_HD[ch][wvf_shift].begin(), min_it, [npe_anode](auto a, auto b) { return a+npe_anode*b; }); } @@ -474,7 +560,7 @@ namespace opdet { } - void DigiPMTSBNDAlg::AddLineNoise(std::vector& wave) + void DigiPMTSBNDAlg::AddLineNoise(std::vector& wave, int ch) { // TODO: after running the profiler I can see that this is where // most cycles are being used. Potentially some improvement could @@ -495,7 +581,66 @@ namespace opdet { } - void DigiPMTSBNDAlg::AddDarkNoise(std::vector& wave) + void DigiPMTSBNDAlg::AddDataNoise(std::vector& wave, int ch) + { + // Get the noise waveform from data + // Get all the noise waveforms that correspond to the channel we are using + // Select a random noise waveform + // Compare the length of the noise and the simulation waveform: + // if the simulation length is larger then choose another random noise waveform to fill up the missing items + + if (!noise_file || noise_file->IsZombie()) { + throw cet::exception("DigiPMTSBNDAlg") << " PMT Noise file could not be opened " << std::endl; + return; + } + + TList *keys = noise_file->GetListOfKeys(); + if (!keys || keys->IsEmpty()) { + throw cet::exception("DigiPMTSBNDAlg") << " PMT Noise file is empty " << std::endl; + return; + } + + std::string opChannelName = "opchannel_" + std::to_string(ch) + "_pmt"; + std::vector keylist; + TIter nextkey(keys); + TKey *key; + while ((key = (TKey*)nextkey())) { + std::string channel_name = key->GetName(); + if (channel_name.find(opChannelName) != std::string::npos) { + keylist.push_back(key); + } + } + + if (keylist.empty()) { + throw cet::exception("DigiPMTSBNDAlg") << " PMT Noise file has no data for channel " << ch << std::endl; + return; + } + + std::vector noise_wform; + int waveBins = wave.size(); + int currentSize = 0; + + // Fill noise vector until it matches the waveform length + while (currentSize < waveBins) { + int noiseWformIdx = static_cast(fEngine->flat() * keylist.size()); + TH1 *noiseHist = (TH1*)keylist[noiseWformIdx]->ReadObj(); + + for (int i = 1; i <= noiseHist->GetNbinsX(); i++) { + noise_wform.push_back(noiseHist->GetBinContent(i)); + currentSize += 1; + if (currentSize >= waveBins) break; + } + delete noiseHist; + } + + // Add noise to the waveform + for (size_t i = 0; i < wave.size(); i++) { + wave[i] += noise_wform[i]; + } + } + + + void DigiPMTSBNDAlg::AddDarkNoise(std::vector& wave, int ch) { double timeBin; // Multiply by 10^9 since fParams.DarkNoiseRate is in Hz (conversion from s to ns) @@ -503,7 +648,7 @@ namespace opdet { double darkNoiseTime = fExponentialGen.fire(mean); while(darkNoiseTime < wave.size()) { timeBin = std::round(darkNoiseTime); - if(timeBin < wave.size()) {AddSPE(fSamplingPeriod*timeBin, wave);} + if(timeBin < wave.size()) {AddSPE(fSamplingPeriod*timeBin, wave, ch);} // Find next time to add dark noise darkNoiseTime += fExponentialGen.fire(mean); } @@ -541,7 +686,6 @@ namespace opdet { return t_min; } - // TODO: this function is not being used anywhere! ~icaza double DigiPMTSBNDAlg::FindMinimumTimeLite( sim::SimPhotonsLite const& litesimphotons, @@ -605,6 +749,8 @@ namespace opdet { fBaseConfig.TTS = config.tts(); fBaseConfig.CableTime = config.cableTime(); fBaseConfig.PMTDataFile = config.pmtDataFile(); + fBaseConfig.UseDataNoise = config.UseDataNoise(); + fBaseConfig.OpDetNoiseFile = config.OpDetNoiseFile(); fBaseConfig.MakeGainFluctuations = config.gainFluctuationsParams.get_if_present(fBaseConfig.GainFluctuationsParams); fBaseConfig.SimulateNonLinearity = config.nonLinearityParams.get_if_present(fBaseConfig.NonLinearityParams); config.hdOpticalWaveformParams.get_if_present(fBaseConfig.HDOpticalWaveformParams); diff --git a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh index ed9190dd4..7a7d7bcd0 100644 --- a/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh +++ b/sbndcode/OpDetSim/DigiPMTSBNDAlg.hh @@ -37,12 +37,19 @@ #include "lardata/DetectorInfoServices/DetectorClocksServiceStandard.h" #include "lardataobj/Simulation/SimPhotons.h" #include "lardata/DetectorInfoServices/LArPropertiesService.h" +#include "larcore/Geometry/WireReadout.h" #include "sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh" #include "sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh" #include "sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms.hh" +#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h" + #include "TFile.h" +#include "TKey.h" +#include "TH1.h" +#include "TH1D.h" namespace opdet { @@ -69,10 +76,14 @@ namespace opdet { double PMTCoatedVISEff_tpc1; //PMT (coated) efficiency for reflected (VIS) light (TPC1) double PMTUncoatedEff_tpc1; //PMT (uncoated) efficiency (TPC1) std::string PMTDataFile; //File containing timing emission structure for TPB, and single PE profile from data - bool PMTSinglePEmodel; //Model for single pe response, false for ideal, true for test bench meas + std::string OpDetNoiseFile; + std::string PMTSinglePEmodel; //Model for single pe response, false for ideal, true for test bench meas bool MakeGainFluctuations; //Fluctuate PMT gain fhicl::ParameterSet GainFluctuationsParams; bool SimulateNonLinearity; //Fluctuate PMT gain + bool PositivePolarity; + bool OscillateAfterPulse; + bool UseDataNoise; fhicl::ParameterSet NonLinearityParams; fhicl::ParameterSet HDOpticalWaveformParams; @@ -124,6 +135,8 @@ namespace opdet { return fParams.PMTBaseline; } + void AddOscillationAfterPulse( std::vector& wave); + private: ConfigurationParameters_t fParams; @@ -138,10 +151,14 @@ namespace opdet { double fPMTUncoatedEff_tpc1; bool fPositivePolarity; int fADCSaturation; + bool fUseDataNoise; + std::string fOpDetNoiseFile; double sigma1; double sigma2; + TFile* noise_file; + const double transitTimeSpread_frac = 2.0 * std::sqrt(2.0 * std::log(2.0)); CLHEP::HepRandomEngine* fEngine; //!< Reference to art-managed random-number engine @@ -151,6 +168,9 @@ namespace opdet { CLHEP::RandExponential fExponentialGen; std::unique_ptr fTimeTPB; // histogram for getting the TPB emission time for coated PMTs + //PMTCalibrationDatabase service + sbndDB::PMTCalibrationDatabase const* fPMTCalibrationDatabaseService; + //PMTFluctuationsAlg std::unique_ptr fPMTGainFluctuationsPtr; //HDWaveforms @@ -159,15 +179,19 @@ namespace opdet { //PMTNonLinearity std::unique_ptr fPMTNonLinearityPtr; - void AddSPE(size_t time, std::vector& wave, double npe = 1); // add single pulse to auxiliary waveform + void AddSPE(size_t time, std::vector& wave, int ch ,double npe = 1); // add single pulse to auxiliary waveform void Pulse1PE(std::vector& wave); double Transittimespread(double fwhm); std::vector fSinglePEWave; // single photon pulse vector - std::vector> fSinglePEWave_HD; // single photon pulse vector + std::vector fAverageDataSER; // single photon pulse vector + double fAverageDataSPEAmplitude; + std::vector>> fSinglePEWave_HD; // single photon pulse vector int pulsesize; //size of 1PE waveform std::unordered_map< raw::Channel_t, std::vector > fFullWaveforms; + geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); + void CreatePDWaveformUncoatedPMT( sim::SimPhotons const& SimPhotons, double t_min, @@ -193,8 +217,9 @@ namespace opdet { std::unordered_map& DirectPhotonsMap, std::unordered_map& ReflectedPhotonsMap); void CreateSaturation(std::vector& wave);//Including saturation effects (dynamic range) - void AddLineNoise(std::vector& wave); //add noise to baseline - void AddDarkNoise(std::vector& wave); //add dark noise + void AddLineNoise(std::vector& wave, int ch); //add noise to baseline + void AddDataNoise(std::vector& wave, int ch); //add noise from data to baseline + void AddDarkNoise(std::vector& wave, int ch); //add dark noise double FindMinimumTime( sim::SimPhotons const&, int ch, @@ -299,9 +324,9 @@ namespace opdet { Comment("PMT (uncoated) detection efficiency (TPC1)") }; - fhicl::Atom PMTsinglePEmodel { + fhicl::Atom PMTsinglePEmodel { Name("PMTSinglePEmodel"), - Comment("Model used for single PE response of PMT. =0 is ideal, =1 is testbench") + Comment("Model used for single PE response of PMT. Accepted strings are \"testbench\", \"data\", and \"ideal\"") }; fhicl::Atom pmtDataFile { @@ -309,6 +334,16 @@ namespace opdet { Comment("File containing timing emission distribution for TPB and single pe pulse from data") }; + fhicl::Atom UseDataNoise { + Name("UseDataNoise"), + Comment("Use data noise from file") + }; + + fhicl::Atom OpDetNoiseFile { + Name("OpDetNoiseFile"), + Comment("Use data noise from this file") + }; + fhicl::OptionalDelegatedParameter gainFluctuationsParams { Name("GainFluctuationsParams"), Comment("Parameters used for SinglePE response fluctuations") @@ -324,6 +359,7 @@ namespace opdet { Comment("Parameters used for high definition waveform") }; + }; //struct Config DigiPMTSBNDAlgMaker(Config const& config); //Constructor diff --git a/sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms_config.fcl b/sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms_config.fcl index 7ee6220b5..5a5e2a23f 100644 --- a/sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms_config.fcl +++ b/sbndcode/OpDetSim/HDWvf/HDOpticalWaveforms_config.fcl @@ -11,7 +11,7 @@ IncludeHDOpticalWaveforms_PMT: { tool_type: "HDOpticalWaveforms" - ResolutionIncrease: 2. + ResolutionIncrease: 1. } END_PROLOG diff --git a/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt b/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt index 115aa504d..e97172251 100644 --- a/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt +++ b/sbndcode/OpDetSim/PMTAlg/CMakeLists.txt @@ -5,6 +5,13 @@ cet_build_plugin(PMTGainFluctuations1Dynode art::tool nurandom::RandomUtils_NuRandomService_service ) +cet_build_plugin(PMTGaussianGainFluctuation art::tool + SOURCE + PMTGaussianGainFluctuation_tool.cc + LIBRARIES + nurandom::RandomUtils_NuRandomService_service + larcore::Geometry_Geometry_service +) cet_build_plugin(PMTNonLinearityTF1 art::tool SOURCE @@ -13,8 +20,16 @@ cet_build_plugin(PMTNonLinearityTF1 art::tool ROOT::Hist ) +cet_build_plugin(PMTNonLinearityTF1ChannelByChannel art::tool + SOURCE + PMTNonLinearityTF1ChannelByChannel_tool.cc + LIBRARIES + art::Framework_Services_Registry + ROOT::Hist + larcore::Geometry_Geometry_service +) install_headers() install_fhicl() install_source() -FILE(GLOB fcl_files *.fcl) +FILE(GLOB fcl_files *.fcl) \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh b/sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh index 6d5a635f9..1c111e593 100644 --- a/sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh +++ b/sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh @@ -20,7 +20,14 @@ public: virtual ~PMTGainFluctuations() noexcept = default; //Returns fluctuated factor for SPR - virtual double GainFluctuation(unsigned int npe, CLHEP::HepRandomEngine* eng) = 0; + virtual double GainFluctuation(unsigned int npe, CLHEP::HepRandomEngine* eng){ + // Default implementation, can be overridden + return npe; + } + virtual double GainFluctuation(int ch, unsigned int npe, CLHEP::HepRandomEngine* eng) { + // Default implementation, can be overridden + return GainFluctuation(npe, eng); + } }; #endif diff --git a/sbndcode/OpDetSim/PMTAlg/PMTGaussianGainFluctuation_tool.cc b/sbndcode/OpDetSim/PMTAlg/PMTGaussianGainFluctuation_tool.cc new file mode 100644 index 000000000..c413cd109 --- /dev/null +++ b/sbndcode/OpDetSim/PMTAlg/PMTGaussianGainFluctuation_tool.cc @@ -0,0 +1,62 @@ +//////////////////////////////////////////////////////////////////////// +// Specific class tool for PMTGainFluctuations +// File: PMTGaussianGainFluctuation.hh +// Base class: PMTGainFluctuations.hh +// Algorithm based on function +// 'multiplicationStageGain(unsigned int i /* = 1 */) const' +// in icaruscode/PMT/Algorithms/PMTsimulationAlg.cxx +//////////////////////////////////////////////////////////////////////// + +#include "fhiclcpp/ParameterSet.h" +#include "art/Utilities/ToolMacros.h" +#include "art/Utilities/make_tool.h" +#include "art/Utilities/ToolConfigTable.h" +#include "nurandom/RandomUtils/NuRandomService.h" +#include "CLHEP/Random/RandGaussQ.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" + +#include + +#include "sbndcode/OpDetSim/PMTAlg/PMTGainFluctuations.hh" +#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h" + +namespace opdet { + class PMTGaussianGainFluctuation; +} + + +class opdet::PMTGaussianGainFluctuation : opdet::PMTGainFluctuations { +public: + + //Configuration parameters + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + }; + + explicit PMTGaussianGainFluctuation(art::ToolConfigTable const& config); + + //Returns fluctuated factor for SPR + double GainFluctuation(int ch, unsigned int npe, CLHEP::HepRandomEngine* eng) override; + +private: + + //PMTCalibrationDatabase service + sbndDB::PMTCalibrationDatabase const* fPMTCalibrationDatabaseService; +}; + + +opdet::PMTGaussianGainFluctuation::PMTGaussianGainFluctuation(art::ToolConfigTable const& config) +{ + fPMTCalibrationDatabaseService = lar::providerFrom(); +} + +double opdet::PMTGaussianGainFluctuation::GainFluctuation(int ch, unsigned int npe, CLHEP::HepRandomEngine* eng){ + double spe_amp = fPMTCalibrationDatabaseService->getSPEAmplitude(ch); + double ChannelGainFluctuation = fPMTCalibrationDatabaseService->getSPEAmplitudeStd(ch); + return CLHEP::RandGaussQ::shoot(eng, npe, std::sqrt(npe)*(ChannelGainFluctuation/spe_amp)); +} + +DEFINE_ART_CLASS_TOOL(opdet::PMTGaussianGainFluctuation) diff --git a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh index 13b4bf04e..5b00f3b56 100644 --- a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh +++ b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh @@ -18,9 +18,19 @@ class opdet::PMTNonLinearity { public: //Constructor virtual ~PMTNonLinearity() noexcept = default; + + //Configure non-linearity tool + virtual void ConfigureNonLinearity() = 0; //Returns rescaled number of PE - virtual double NObservedPE(size_t bin, std::vector & pe_vector) = 0; + virtual double NObservedPE(size_t bin, std::vector & pe_vector){ + return pe_vector[bin]; + } + virtual double NObservedPE(int opch, size_t bin, std::vector & pe_vector){ + // Default implementation calls the non-channel-dependent version + return NObservedPE(bin, pe_vector); + } + }; #endif diff --git a/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc new file mode 100644 index 000000000..c7e217b91 --- /dev/null +++ b/sbndcode/OpDetSim/PMTAlg/PMTNonLinearityTF1ChannelByChannel_tool.cc @@ -0,0 +1,125 @@ +//////////////////////////////////////////////////////////////////////// +// Specific class tool for PMTNonLinearity +// File: PMTNonLinearityTF1_tool.hh +// Base class: PMTNonLinearity.hh +//////////////////////////////////////////////////////////////////////// + +#include "fhiclcpp/ParameterSet.h" +#include "art/Utilities/ToolMacros.h" +#include "art/Utilities/ToolConfigTable.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" + +#include "larcore/Geometry/WireReadout.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" + +#include +#include +#include +#include "TF1.h" + +#include "sbndcode/OpDetSim/PMTAlg/PMTNonLinearity.hh" +#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h" + +namespace opdet { + class PMTNonLinearityTF1ChannelByChannel; +} + + +class opdet::PMTNonLinearityTF1ChannelByChannel : opdet::PMTNonLinearity { +public: + + //Configuration parameters + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + + fhicl::Atom attenuationForm { + Name("AttenuationForm"), + Comment("Non linearity functional form") + }; + + fhicl::Atom attenuationPreTime { + Name("AttenuationPreTime"), + Comment("For the non-linear attenuation, consider photons arriving given by this time window. In nanoseconds.") + }; + + fhicl::Sequence nonLinearRange { + Name("NonLinearRange"), + Comment("Non linear range. Assume linear response/completely saturated response out of this range. In #PE.") + }; + + }; + + explicit PMTNonLinearityTF1ChannelByChannel(art::ToolConfigTable const& config); + + //Returns rescaled #pe after non linearity + double NObservedPE(int opch, size_t bin, std::vector & pe_vector) override; + void ConfigureNonLinearity() override; + +private: + //Configuration parameters + std::string fAttenuationForm; + unsigned int fAttenuationPreTime; + std::vector fNonLinearRange; + geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); + + //TF1 for non linearity function + std::map> fNonLinearTF1Map; + + // Vector to store non linearity attenuation values + std::vector> fPEAttenuation_V; + std::vector fPEAttenuation; + std::vector fPESaturationValue_V; + int fPESaturationValue; + + //PMTCalibrationDatabase service + sbndDB::PMTCalibrationDatabase const* fPMTCalibrationDatabaseService; +}; + + +opdet::PMTNonLinearityTF1ChannelByChannel::PMTNonLinearityTF1ChannelByChannel(art::ToolConfigTable const& config) + : fAttenuationForm { config().attenuationForm() } + , fAttenuationPreTime { config().attenuationPreTime() } + , fNonLinearRange { config().nonLinearRange() } +{ + fPMTCalibrationDatabaseService = lar::providerFrom(); +} + + +void opdet::PMTNonLinearityTF1ChannelByChannel::ConfigureNonLinearity(){ + + //Initialize TF1 for each channel with channel-dependent parameters + for(size_t opch=0; opch(("NonLinearTF1_"+std::to_string(opch)).c_str(), fAttenuationForm.c_str()); + fNonLinearTF1Map[opch]->SetParameter(0, fPMTCalibrationDatabaseService->getNonLineatiryPESat(opch)); + fNonLinearTF1Map[opch]->SetParameter(1, fPMTCalibrationDatabaseService->getNonLineatiryAlpha(opch)); + } + + // Initialize attenuation vector + fPEAttenuation_V.resize(fWireReadout.NOpChannels()); + for(size_t opch=0; opchEval(pe)/pe; + } + } + + fPESaturationValue_V.resize(fWireReadout.NOpChannels(), 1); + for(size_t opch=0; opchEval(fNonLinearRange[1])); + } +} + + +double opdet::PMTNonLinearityTF1ChannelByChannel::NObservedPE(int opch, size_t bin, std::vector & pe_vector){ + size_t start_bin = bin-fAttenuationPreTime; + if(fAttenuationPreTime<0) start_bin=0; + unsigned int npe_acc = std::accumulate(pe_vector.begin()+start_bin,pe_vector.begin()+bin+1, 0); + + if(npe_acc & pe_vector) override; + void ConfigureNonLinearity() override; private: //Configuration parameters @@ -66,7 +67,7 @@ class opdet::PMTNonLinearityTF1 : opdet::PMTNonLinearity { std::vector fNonLinearRange; //TF1 for non linearity function - TF1 *fNonLinearTF1; + std::unique_ptr fNonLinearTF1; // Vector to store non linearity attenuation values std::vector fPEAttenuation_V; @@ -79,8 +80,10 @@ opdet::PMTNonLinearityTF1::PMTNonLinearityTF1(art::ToolConfigTable const , fAttenuationFormParams { config().attenuationFormParams() } , fAttenuationPreTime { config().attenuationPreTime() } , fNonLinearRange { config().nonLinearRange() } -{ - fNonLinearTF1 = new TF1("NonLinearTF1", fAttenuationForm.c_str()); +{} + +void opdet::PMTNonLinearityTF1::ConfigureNonLinearity(){ + fNonLinearTF1 = std::make_unique("NonLinearTF1", fAttenuationForm.c_str()); for(size_t k=0; kSetParameter(k, fAttenuationFormParams[k]); } diff --git a/sbndcode/OpDetSim/PMTAlg/pmtgainfluctuations_config.fcl b/sbndcode/OpDetSim/PMTAlg/pmtgainfluctuations_config.fcl index 5b611b525..59e59b535 100644 --- a/sbndcode/OpDetSim/PMTAlg/pmtgainfluctuations_config.fcl +++ b/sbndcode/OpDetSim/PMTAlg/pmtgainfluctuations_config.fcl @@ -18,4 +18,9 @@ FirstDynodeGainFluctuations: Gain: 1e7 # total typical PMT gain } +GaussianGainFluctuations: +{ + tool_type: "PMTGaussianGainFluctuation" +} + END_PROLOG diff --git a/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl b/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl index 82c6e3c53..5fdeca0ef 100644 --- a/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl +++ b/sbndcode/OpDetSim/PMTAlg/pmtnonlinearity_config.fcl @@ -12,4 +12,14 @@ PMTNonLinearityTF1: NonLinearRange: [100, 701] } +PMTNonLinearityTF1ChannelByChannel: +{ + tool_type: "PMTNonLinearityTF1ChannelByChannel" + #AttenuationForm: "sqrt( (sqrt(1+8*(x/[0])**[1])-1)/( 4*(x/[0])**[1] ) )" + #AttenuationFormParams: [78.15, 1.96] + AttenuationForm: "x/sqrt(1+(x/[0])**[1])" + AttenuationPreTime: 4 #ns + NonLinearRange: [100, 7000] +} + END_PROLOG diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/CMakeLists.txt b/sbndcode/OpDetSim/PMTPulseOscillations/CMakeLists.txt new file mode 100644 index 000000000..2eb94cd3f --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/CMakeLists.txt @@ -0,0 +1,28 @@ +set( + MODULE_LIBRARIES + sbndcode_OpDetSim + larcore::Geometry_Geometry_service + lardataobj::Simulation + lardata::Utilities + lardataobj::RawData + sbndcode_Utilities_SignalShapingServiceSBND_service + nurandom::RandomUtils_NuRandomService_service + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Optional_RandomNumberGenerator_service + canvas::canvas + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + cetlib::cetlib + cetlib_except::cetlib_except + CLHEP::CLHEP + ${ROOT_BASIC_LIB_LIST} +) + +cet_build_plugin(PMTPulseOscillation art::module SOURCE PMTPulseOscillation_module.cc LIBRARIES ${MODULE_LIBRARIES}) + +install_headers() +install_fhicl() +install_source() +add_subdirectory(job) +##FILE(GLOB fcl_files *.fcl) \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc b/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc new file mode 100644 index 000000000..3f514fb17 --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/PMTPulseOscillation_module.cc @@ -0,0 +1,145 @@ +//////////////////////////////////////////////////////////////////////// +// Class: SBNDOpDeconvolution +// Plugin Type: producer (art v3_06_03) +// File: SBNDOpDeconvolution_module.cc +// +// Generated at Tue Jul 13 06:29:02 2021 by Francisco Nicolas-Arnaldos using cetskelgen +// from cetlib version v3_11_01. +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "art/Utilities/make_tool.h" +#include "art_root_io/TFileService.h" + +#include "CLHEP/Random/RandFlat.h" +#include "nurandom/RandomUtils/NuRandomService.h" + +#include + +#include "lardataobj/RawData/OpDetWaveform.h" + +#include "sbndcode/OpDetSim/sbndPDMapAlg.hh" +#include "TF1.h" +#include "TH1F.h" + +namespace opdet { + class PMTPulseOscillation; +} + + +class opdet::PMTPulseOscillation : public art::EDProducer { + public: + explicit PMTPulseOscillation(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + PMTPulseOscillation(PMTPulseOscillation const&) = delete; + PMTPulseOscillation(PMTPulseOscillation&&) = delete; + PMTPulseOscillation& operator=(PMTPulseOscillation const&) = delete; + PMTPulseOscillation& operator=(PMTPulseOscillation&&) = delete; + + // Required functions. + void produce(art::Event& e) override; + void CreateOscillatedWaveform(const raw::OpDetWaveform& oldWaveform , raw::OpDetWaveform& newWaveform ); + + + private: + + // Declare member data here. + std::string fInputLabel; + std::vector fPDTypes; + std::vector fElectronics; + std::string fOscillationFunction; + double fOscillationFrequency; + double fOscillationDampingConstant; + double fSamplingPeriod; + double fPulseOscillationThreshold; + int fOscillationOffset; + //PDS map + opdet::sbndPDMapAlg pdsmap; + std::vector fOscillationFunctionParams; + TF1 *fOscillationTF1; + CLHEP::HepRandomEngine& fEngine; +}; + + +opdet::PMTPulseOscillation::PMTPulseOscillation(fhicl::ParameterSet const& p) + : EDProducer{p}, + fEngine(art::ServiceHandle{}->registerAndSeedEngine( + createEngine(0, "HepJamesRandom", "pulseosc"), "HepJamesRandom", "pulseosc", p, "Seed")) + // More initializers here. +{ + // Call appropriate produces<>() functions here. + // Call appropriate consumes<>() for any products to be retrieved by this module. + fInputLabel = p.get< std::string >("InputLabel"); + fPDTypes = p.get< std::vector >("PDTypes"); + fElectronics = p.get< std::vector >("Electronics"); + fOscillationFunction = p.get("OscillationFunction"); + fSamplingPeriod = p.get("SamplingPeriod"); + fOscillationFrequency = p.get("OscillationFrequency"); + fOscillationDampingConstant = p.get("OscillationDampingConstant"); + fOscillationFunctionParams = p.get>("OscillationFunctionParams"); + fPulseOscillationThreshold = p.get("PulseOscillationThreshold"); + fOscillationOffset = p.get("OscillationOffset"); + fOscillationTF1 = new TF1("OscillationTemplate", fOscillationFunction.c_str()); + for(size_t i=0; iSetParameter(i, fOscillationFunctionParams[i]); + } + + produces< std::vector< raw::OpDetWaveform > >(); +} + +void opdet::PMTPulseOscillation::produce(art::Event& e) +{ + //Load the waveforms + art::Handle< std::vector< raw::OpDetWaveform > > wfHandle; + e.getByLabel(fInputLabel, wfHandle); + if (!wfHandle.isValid()) { + std::cout << " not finding the handles " << std::endl; + mf::LogError("PMTPulseOscillation")<<"Input waveforms with input label "< > OscillatedWf_VecPtr(std::make_unique< std::vector< raw::OpDetWaveform > > ()); + for(auto const& wf : *wfHandle){ + if(std::find(fPDTypes.begin(), fPDTypes.end(), pdsmap.pdType(wf.ChannelNumber()) ) != fPDTypes.end() ){ + raw::OpDetWaveform new_wf = wf; + CreateOscillatedWaveform(wf, new_wf); + OscillatedWf_VecPtr->push_back(new_wf); + } + else OscillatedWf_VecPtr->push_back(wf); + } + + e.put( std::move(OscillatedWf_VecPtr) ); +} + + +void opdet::PMTPulseOscillation::CreateOscillatedWaveform(const raw::OpDetWaveform& oldWaveform, raw::OpDetWaveform& newWaveform) +{ + //Find the idx of the maximum (or minimum) point of the waveform + int starting_ttick = std::distance(oldWaveform.begin(), std::min_element(oldWaveform.begin(), oldWaveform.end())); + double baseline = accumulate( oldWaveform.begin(), oldWaveform.end(), 0.0) / oldWaveform.size(); + double pulse_max = abs(oldWaveform[starting_ttick]-baseline); + if(pulse_maxEval(pulse_max); + for(size_t i=starting_ttick; i(i)-static_cast(starting_ttick))*fSamplingPeriod; + double envelope = std::exp(-delta_time / fOscillationDampingConstant); + double phase = 2.0 * M_PI * fOscillationFrequency * delta_time; + double new_value = static_cast(newWaveform.Waveform()[i+fOscillationOffset]) + oscillation_amplitude * envelope * std::cos(phase+M_PI); + // Add some dither to avoid quantization effects + double dither = CLHEP::RandFlat::shoot(&fEngine,1.0) - 0.5; + newWaveform.Waveform()[i+fOscillationOffset] = std::round(new_value + dither); + if(envelope < 0.05) return; + } +} + +DEFINE_ART_MODULE(opdet::PMTPulseOscillation) \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/job/CMakeLists.txt b/sbndcode/OpDetSim/PMTPulseOscillations/job/CMakeLists.txt new file mode 100644 index 000000000..eca2a52d6 --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/job/CMakeLists.txt @@ -0,0 +1 @@ +install_fhicl() \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl b/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl new file mode 100644 index 000000000..4e727f6dd --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/job/run_pmtpulseoscillation.fcl @@ -0,0 +1,62 @@ +#include "services_sbnd.fcl" +#include "simulationservices_sbnd.fcl" +#include "messages_sbnd.fcl" +#include "sam_sbnd.fcl" +#include "sbnd_pmtpulseoscillation_config.fcl" + + +process_name: pmtpulseoscillation + +services: +{ + # Load the service that manages root files for histograms. + TFileService: { fileName: "hist.root" } + RandomNumberGenerator: {} # required by fuzzyCluster + message: @local::sbnd_message_services_prod # from messages_sbnd.fcl + NuRandomService: @local::sbnd_seedservice +} + +#source is now a root file +source: +{ + module_type: RootInput + maxEvents: -1 # Number of events to create +} + +outputs: +{ + out1: + { + module_type: RootOutput + fileName: "%ifb_%p.root" + dataTier: "reconstructed" + compressionLevel: 1 + } +} + +physics: +{ + + producers: + { + # random number saver + rns: { module_type: "RandomNumberSaver" } + pmtpulseoscillation: @local::sbnd_pmtpulseoscillation + } + + + # define the producer and filter modules for this path, order matters, + # filters reject all following items. see lines starting physics.producers below + reco: [rns, pmtpulseoscillation ] + + stream1: [out1] + + # trigger_paths is a keyword and contains the paths that modify the art::event, + # ie filters and producers + trigger_paths: [reco] + + + # end_paths is a keyword and contains the paths that do not modify the art::Event, + # ie analyzers and output streams. these all run simultaneously + end_paths: [stream1] +} \ No newline at end of file diff --git a/sbndcode/OpDetSim/PMTPulseOscillations/job/sbnd_pmtpulseoscillation_config.fcl b/sbndcode/OpDetSim/PMTPulseOscillations/job/sbnd_pmtpulseoscillation_config.fcl new file mode 100644 index 000000000..143a9762b --- /dev/null +++ b/sbndcode/OpDetSim/PMTPulseOscillations/job/sbnd_pmtpulseoscillation_config.fcl @@ -0,0 +1,19 @@ +BEGIN_PROLOG + +sbnd_pmtpulseoscillation: +{ + # Parameters for ideal SER simulation + module_type: "PMTPulseOscillation" + InputLabel: "opdaq" + PDTypes: ["pmt_coated", "pmt_uncoated"] + Electronics: [] + OscillationFrequency: 0.000178571 + OscillationDampingConstant: 50000 + SamplingPeriod: 2 + OscillationFunction: "[0]*exp([1] * x)" + OscillationFunctionParams: [0.20292146, 0.00079746] + PulseOscillationThreshold: 70 + OscillationOffset: 0 +} + +END_PROLOG \ No newline at end of file diff --git a/sbndcode/OpDetSim/digi_pmt_sbnd.fcl b/sbndcode/OpDetSim/digi_pmt_sbnd.fcl index a5d6ca7d2..801220e7a 100644 --- a/sbndcode/OpDetSim/digi_pmt_sbnd.fcl +++ b/sbndcode/OpDetSim/digi_pmt_sbnd.fcl @@ -15,8 +15,8 @@ sbnd_digipmt_alg: PMTChargeToADC: -25.97 #charge to adc factor # Parameters for test bench SER simulation - PMTSinglePEmodel: true #false for ideal PMT response, true for test bench measured response - PMTDataFile: "OpDetSim/digi_pmt_sbnd_v2int0.root" # located in sbnd_data + PMTSinglePEmodel: "data" + PMTDataFile: "OpDetSim/digi_pmt_sbnd_v2int0.root" # located in sbnd_data # Time delays TTS: 2.4 #Transit Time Spread in ns @@ -27,25 +27,29 @@ sbnd_digipmt_alg: PMTBaseline: 14745 #in ADC PMTBaselineRMS: 2.6 #in ADC + #Sample noise from data + UseDataNoise: true + OpDetNoiseFile: "OpDetReco/PMTNoiseTemplate.root" + # Dark counts PMTDarkNoiseRate: 1000.0 #in Hz # Detection efficiencies - PMTCoatedVUVEff_tpc0: 0.035 #PMT coated detection efficiency for direct (VUV) light - PMTCoatedVISEff_tpc0: 0.03882 #PMT coated detection efficiency for reflected (VIS) light - PMTUncoatedEff_tpc0: 0.03758 #PMT uncoated detection efficiency + PMTCoatedVUVEff_tpc0: 0.0315 #PMT coated detection efficiency for direct (VUV) light + PMTCoatedVISEff_tpc0: 0.03493 #PMT coated detection efficiency for reflected (VIS) light + PMTUncoatedEff_tpc0: 0.03382 #PMT uncoated detection efficiency - PMTCoatedVUVEff_tpc1: 0.03 #PMT coated detection efficiency for direct (VUV) light - PMTCoatedVISEff_tpc1: 0.03882 #PMT coated detection efficiency for reflected (VIS) light - PMTUncoatedEff_tpc1: 0.03758 #PMT uncoated detection efficiency + PMTCoatedVUVEff_tpc1: 0.027 #PMT coated detection efficiency for direct (VUV) light + PMTCoatedVISEff_tpc1: 0.03493 #PMT coated detection efficiency for reflected (VIS) light + PMTUncoatedEff_tpc1: 0.03382 #PMT uncoated detection efficiency # Simulate gain fluctuations # (comment-out to skip gain fluctuations simulation) - GainFluctuationsParams: @local::FirstDynodeGainFluctuations + GainFluctuationsParams: @local::GaussianGainFluctuations # Simulate PMT non linear response # (comment-out to skip non linearities simulation) - NonLinearityParams: @local::PMTNonLinearityTF1 + NonLinearityParams: @local::PMTNonLinearityTF1ChannelByChannel HDOpticalWaveformParamsPMT: @local::IncludeHDOpticalWaveforms_PMT }