From 8fcebb4c49e0c3d3eb59f212e16dc6b3bfdc9b2a Mon Sep 17 00:00:00 2001 From: Nikki Pallat Date: Mon, 22 Sep 2025 10:46:42 -0500 Subject: [PATCH 1/7] Apply TriggerWorkNikki changes on top of develop --- .../LArSoftConfigurations/services_sbnd.fcl | 2 + sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl | 47 ++ sbndcode/OpDetSim/BeamRatesCalib_module.cc | 793 ++++++++++++++++++ sbndcode/OpDetSim/TriggerEmulationService.cc | 191 +++++ sbndcode/OpDetSim/TriggerEmulationService.h | 63 ++ .../TriggerEmulationService_service.cc | 4 + .../OpDetSim/opDetDigitizerSBND_module.cc | 416 ++++++++- sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl | 7 + sbndcode/OpDetSim/run_BeamRatesCalib.fcl | 37 + .../OpDetSim/triggeremulationservice_sbnd.fcl | 15 + 10 files changed, 1543 insertions(+), 32 deletions(-) create mode 100644 sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl create mode 100755 sbndcode/OpDetSim/BeamRatesCalib_module.cc create mode 100644 sbndcode/OpDetSim/TriggerEmulationService.cc create mode 100644 sbndcode/OpDetSim/TriggerEmulationService.h create mode 100644 sbndcode/OpDetSim/TriggerEmulationService_service.cc mode change 100644 => 100755 sbndcode/OpDetSim/opDetDigitizerSBND_module.cc mode change 100644 => 100755 sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl create mode 100644 sbndcode/OpDetSim/run_BeamRatesCalib.fcl create mode 100644 sbndcode/OpDetSim/triggeremulationservice_sbnd.fcl diff --git a/sbndcode/LArSoftConfigurations/services_sbnd.fcl b/sbndcode/LArSoftConfigurations/services_sbnd.fcl index 6ba1cfe07..a14322b2f 100644 --- a/sbndcode/LArSoftConfigurations/services_sbnd.fcl +++ b/sbndcode/LArSoftConfigurations/services_sbnd.fcl @@ -46,6 +46,7 @@ #include "sam_sbnd.fcl" #include "spacecharge_sbnd.fcl" #include "pmtcalibrationdatabase_sbnd.fcl" +#include "triggeremulationservice_sbnd.fcl" BEGIN_PROLOG @@ -90,6 +91,7 @@ sbnd_services: DetPedestalService: @local::sbnd_detpedestalservice # from database_sbnd.fcl SpaceCharge: @local::sbnd_spacecharge IPMTCalibrationDatabaseService: @local::sbnd_pmtcalibrationdatabaseservice + TriggerEmulationService: @local::sbnd_triggeremulation_service } # sbnd_services diff --git a/sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl b/sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl new file mode 100644 index 000000000..519ad240b --- /dev/null +++ b/sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl @@ -0,0 +1,47 @@ +BEGIN_PROLOG + +BeamRatesCalib: +{ + module_type: "BeamRatesCalib" + ClusterModuleLabel: "linecluster" + OpDetsToPlot: ["pmt_coated", "pmt_uncoated"] + DaphneFrequency: 62.5 #in MHz. Frequency of the Daphne Readouts + # Width of CAEN trigger pulse in CAEN clock ticks (8ns) + MonWidth: 12 #in 8ns bins + # Stepped over per-PMT MON thresholds for contribution to trigger pulse + MonStart: 24 + MonStop: 64 + MonStep: 1 + # Plotting option for histograms generated + MTCAStart: 1 + MTCAStop: 63 + MTCAStep: 1 + # Save some example trigger pulses for each flash (TH1D similar to wvfana) + SaveAllMon: true + # The MON threshold used in the saved TH1D examples + FCLthreshold: 25 + # PMT index corresponding to the opening of the beam acceptance BeamWindowEnd + # Note that BeamWindowStart should not account for MTC/A wire+processing delays + # Those are handled by CAENOffset + BeamWindowStart: 1000 + BeamWindowEnd: 2500 + CAENOffset: -150 #Based on peak index for flash triggered readouts. Accounts for cable + processing delays + # Sample length of a normal PMT readout + NominalWaveformSize: 5000 + NominalGoodStartTime: 10 #ns. extra offset from TDC + # Flags for additional outputs for MSUM digitization + software trigger + CheckHardwareTriggers: true + CheckSoftTrig: true + # Number of CAEN boards on for this run + TotalCAENBoards: 10 + # Assorted labels to get data products out of Art files + SoftTrigLabel: "pmtmetricproducer" + InputModule: "pmtdecoder" + InputProcess: "DECODE" + InputInstance: "PMTChannels" + PTBLabel: "ptbdecoder::DECODE" + TimingInstanceName: "pmtdecoder:TimingChannels:DECODE" + +} + +END_PROLOG diff --git a/sbndcode/OpDetSim/BeamRatesCalib_module.cc b/sbndcode/OpDetSim/BeamRatesCalib_module.cc new file mode 100755 index 000000000..a6ee5bc75 --- /dev/null +++ b/sbndcode/OpDetSim/BeamRatesCalib_module.cc @@ -0,0 +1,793 @@ +//////////////////////////////////////////////////////////////////////// +// Class: BeamRatesCalib +// Module Type: analyzer +// File: BeamRatesCalib.cc +// +// Analyzer taking in every-beam-spill data and determines rate of flash coincidences with scanned MON thresholds and MTC/A threshold +// +// Authors: J. McLaughlin; Based off wvfAna_module.cc +//////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include + +#include "art/Framework/Core/EDAnalyzer.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 "art_root_io/TFileService.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "canvas/Utilities/Exception.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "larcore/Geometry/Geometry.h" + +#include "lardata/DetectorInfoServices/DetectorClocksService.h" +#include "lardata/DetectorInfoServices/DetectorClocksServiceStandard.h" +#include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "lardata/DetectorInfoServices/LArPropertiesService.h" +#include "lardataobj/RawData/OpDetWaveform.h" +#include "sbndcode/Utilities/SignalShapingServiceSBND.h" +#include "lardataobj/Simulation/sim.h" +#include "lardataobj/Simulation/SimChannel.h" +#include "lardataobj/Simulation/SimPhotons.h" +#include "sbnobj/SBND/Trigger/pmtTrigger.hh" +#include "sbndaq-artdaq-core/Obj/SBND/pmtSoftwareTrigger.hh" +#include "sbndcode/Decoders/PTB/sbndptb.h" + +#include "TH1D.h" +#include "TH2D.h" +#include "TFile.h" +#include "TTree.h" +#include "TGraph.h" + +#include "sbndcode/OpDetSim/sbndPDMapAlg.hh" + +#include "TriggerEmulationService.h" + +namespace opdet { //OpDet means optical detector + + class BeamRateCalib; + + class BeamRateCalib : public art::EDAnalyzer { + public: + explicit BeamRateCalib(fhicl::ParameterSet const & p); + // The destructor generated by the compiler is fine for classes + // without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + //idk why things are doubled here + BeamRateCalib(BeamRateCalib const &) = delete; + BeamRateCalib(BeamRateCalib &&) = delete; + BeamRateCalib & operator = (BeamRateCalib const &) = delete; + BeamRateCalib & operator = (BeamRateCalib &&) = delete; + + //Standard ART loop fucntions + // Required functions. + void analyze(art::Event const & e) override; + //Selected optional functions + void beginJob() override; + void endJob() override; + + void ResetTree(); + void SaveChannelWaveforms(const raw::OpDetWaveform &wvf, int EventCounter, int FlashCounter); + void SaveMonWaveforms(art::Handle< std::vector< raw::OpDetWaveform > > &waveHandle, int MonThreshold, int EventCounter); + void analyzeTrigger(art::Handle< std::vector< raw::OpDetWaveform > > &waveHandle); + void Validate_MonSim_Outputs(art::Event const & e); + double GetWaveformMedian(raw::OpDetWaveform wvf); + opdet::sbndPDMapAlg pdMap; //map for photon detector types + private: + //Data members for analysis/output + //Start with fcl params + int fMonWidth; //Value*10ns MON pulse width when treshold is crossed + //TTree *fWaveformTree; + // Declare member data here. + std::string fInputModuleName; + std::string fInputProcessName; + std::string fInputInstanceName; + std::vector fOpDetsToPlot; //keep for now but probably don't need + std::string opdetType; + std::string fTimingInstanceName; + std::string opdetElectronics; + TH2D* hist_PropTriggers; //Initializes as nullptr + TH2D* hist_PropTriggers_OffBeam; //Initializes as nullptr + TH2D* hist_PropTriggers_NoRestriction; + TH1D* hist_TriggerThreshold; //Initializes as nullptr + //Helpful members to have for the analsis + std::vector MONThresholds; //ADC counts to output a MON pulse + std::vector MTCA_Thresholds; //Number of pairs needed to trigger MTCA output + std::vector hist_TopHatPlots; + TH1D* hist_MSUM_DutyCycle; + TH1D* hist_MSUM_Energy; + TH1D* hist_MSUM_Trigger_Eff; + int EventCounter; + int PMTReadouts; + int PMTPerBoard; + int fChNumber; + int fEvNumber; + int fMonStart; + int fMonStop; + int fMonStep; + int fMTCAStart; + int fMTCAStop; + int fMTCAStep; + int fFCLthreshold; + int fBeamWindowStart; + int fBeamWindowEnd; + bool fSaveAllMON; + bool fCheckTriggers; + bool fMakeTree; + int hist_id; + bool fCheckSoftTrig; + int fNominalGoodStart; + std::string fPTBLabel; + std::string fSoftTrigLabel; + int fTotalCAENBoards; + int fNominalWaveformSize; + bool fTriggerProdCheck; + int fCAENOffset; + int fEtrigOffset; + ULong64_t GateOpenTime; + ULong64_t ETrigTime; + //Event Tree saving infomation + // --- TTrees + TTree* evtTree; + // --- Event information --- + int tree_event; // event number + int tree_run; // run number + int tree_subrun; // subrun number + unsigned int tree_timestamp; // unix time of event + // -- Hardware Trigger Information + std::vector tree_OnBeamMonPeaks; + std::vector tree_OffBeamMonPeaks; + std::vector tree_FullBeamMonPeaks; + std::vector tree_OnBeamMonPeakSample; + std::vector tree_FullBeamMonPeakSample; + int tree_MonStart; + int tree_MonEnd; + int tree_MonStep; + //Software trigger information + int tree_nAboveThreshold; + int tree_trig_ts; // relative to beam window start, in ns + double tree_promptPE; // total PE for the 100 ns following trigger time + double tree_prelimPE; // total PE for the 1 us before trigger time + double tree_peakPE; + double tree_peaktime; + std::vector tree_HLTs; + std::vector tree_HLTTimes; + std::vector tree_LLTs; + std::vector tree_LLTTimes; + + + }; + + BeamRateCalib::BeamRateCalib(fhicl::ParameterSet const & p) + : + EDAnalyzer(p), + hist_PropTriggers(nullptr), + hist_PropTriggers_OffBeam(nullptr), + hist_PropTriggers_NoRestriction(nullptr), + hist_TriggerThreshold(nullptr), + hist_MSUM_DutyCycle(nullptr), + hist_MSUM_Energy(nullptr), + hist_MSUM_Trigger_Eff(nullptr) + // More initializers here. + { + //Read in assorted fcl parameters + GateOpenTime=0; + fInputModuleName = p.get< std::string >("InputModule" ); + fInputProcessName = p.get< std::string >("InputProcess" ); + fInputInstanceName = p.get< std::string >("InputInstance" ); + fTimingInstanceName = p.get("TimingInstanceName"); + fOpDetsToPlot = p.get >("OpDetsToPlot"); + fNominalGoodStart = p.get("NominalGoodStartTime", 2000); + fMonWidth = p.get("MonWidth"); + fMonStart = p.get("MonStart", 20); + fMonStop = p.get("MonStop", 250); + fMonStep = p.get("MonStep", 20); + fMTCAStart = p.get("MTCAStart", 3); + fMTCAStop = p.get("MTCAStop", 64); + fMTCAStep = p.get("MTCAStep", 1); + //PMTPerBoard=15; + PMTPerBoard=16; + fBeamWindowStart = p.get("BeamWindowStart", 828+680); + fBeamWindowEnd = p.get("BeamWindowEnd", 1688+680); + fSaveAllMON = p.get("SaveAllMon", false); + fCheckTriggers = p.get("CheckHardwareTriggers", false); //Needs MTCA LLT to be digitized which is unusual (run 15670) + fFCLthreshold = p.get("FCLthreshold", 15); + fOpDetsToPlot = p.get >("OpDetsToPlot"); + fCheckSoftTrig = p.get("CheckSoftTrig", false); + fSoftTrigLabel = p.get("SoftTrigLabel", "pmtmetricproducer:"); + //fTotalCAENBoards = p.get("TotalCAENBoards", 8); + fTotalCAENBoards = p.get("TotalCAENBoards", 10); + fPTBLabel = p.get< std::string >("PTBLabel", "ptbdecoder::DECODE"); + fTriggerProdCheck = p.get("TriggerProdCheck", false); + fNominalWaveformSize = p.get("NominalWaveformSize", 5000); + fMakeTree = p.get("MakeTree", true); + fCAENOffset = p.get("CAENOffset", 287); + fChNumber = 0; + fEvNumber = 0; + fEtrigOffset = p.get("EtrigOffset", 257); + hist_id=0; + } + + void BeamRateCalib::ResetTree() + { + tree_event = -1; + tree_run = -1; + tree_subrun = -1; + tree_timestamp = -9999; + // -- Hardware Trigger Information + tree_OnBeamMonPeaks.clear(); + tree_OffBeamMonPeaks.clear(); + tree_FullBeamMonPeaks.clear(); + tree_OnBeamMonPeakSample.clear(); + tree_FullBeamMonPeakSample.clear(); + int TotalMonBins = (fMonStop-fMonStart)/fMonStep; + tree_OnBeamMonPeaks.resize(TotalMonBins); + tree_OffBeamMonPeaks.resize(TotalMonBins); + tree_FullBeamMonPeaks.resize(TotalMonBins); + tree_OnBeamMonPeakSample.resize(TotalMonBins); + tree_FullBeamMonPeakSample.resize(TotalMonBins); + tree_HLTs.clear(); + tree_HLTs.resize(200); + tree_HLTTimes.clear(); + tree_HLTTimes.resize(200); + tree_LLTs.clear(); + tree_LLTs.resize(1300); + tree_LLTTimes.clear(); + tree_LLTTimes.resize(1300); + for(int i=0; i tfs; + hist_PropTriggers = tfs->make("hist_PropTriggers", "Proportion of Beam Spill Making Light Trigger", + int(double(MonEnd-MonStart)/MonStep)+1, MonStart-MonStep/2., MonEnd+MonStep/2., int(double(MTCAEnd-MTCAStart)/MTCAStep)+1, MTCAStart-MTCAStep/2.0, MTCAEnd+MTCAStep/2.0); //Give it a real address + bins + hist_PropTriggers_NoRestriction = tfs->make("hist_PropTriggers_LLT", "Proportion of Beam Spills with light LLT", + int(double(MonEnd-MonStart)/MonStep)+1, MonStart-MonStep/2., MonEnd+MonStep/2., int(double(MTCAEnd-MTCAStart)/MTCAStep)+1, MTCAStart-MTCAStep/2.0, MTCAEnd+MTCAStep/2.0); + hist_PropTriggers_OffBeam = tfs->make("hist_PropTriggers_OffBeam", "Proportion of Shifted Beam Spill Making Light Trigger", + int(double(MonEnd-MonStart)/MonStep)+1, MonStart-MonStep/2., MonEnd+MonStep/2., int(double(MTCAEnd-MTCAStart)/MTCAStep)+1, MTCAStart-MTCAStep/2.0, MTCAEnd+MTCAStep/2.0); //Give it a real address + bins + for(int i=0; (MonStart+i*MonStep)<=MonEnd; i++ ) + { + for(int j=0; MTCAStart+j*MTCAStep<=MTCAEnd; j++) + { + std::string histName = "hist_TopHat_Mon_" + std::to_string(MonStart+i*MonStep) + "_MTCA_" + std::to_string(MTCAStart+j*MTCAStep); + std::string histTitle = "hist_TopHat_Mon:" + std::to_string(MonStart+i*MonStep) + "_MTCA:" + std::to_string(MTCAStart+j*MTCAStep); + auto tempHist = tfs->make(histName.c_str(), histTitle.c_str(), 5000, 0, 4999); + hist_TopHatPlots.push_back(tempHist); + int Index = j + i*int((MTCAEnd-MTCAStart)/MTCAStep); + hist_TopHatPlots[Index]->GetXaxis()->SetTitle("Virtual Flash Trigger Index"); + } + } + if(fCheckTriggers) + { + //Make histgoram of MaxPassed for events with a hardware trigger. See if we can pull out the threshold from it + std::string histname = "PMT_Passed_Per_Trigger"; + hist_TriggerThreshold = tfs->make(histname.c_str(), histname.c_str(), 121, -0.5, 120.5); + std::string histname2 = "MSUM_Dutycyle"; + hist_MSUM_DutyCycle = tfs->make(histname2.c_str(),histname2.c_str(), 201, -0.25, 100.25 ); + std::string histname3 = "MSUM_Energy_Joulse"; + hist_MSUM_Energy = tfs->make(histname3.c_str(),histname3.c_str(), 100, 0, 0.000005 ); + std::string histname4 = "MSUM_Trigger_Eff_"; + histname4 = histname4 + std::to_string(fFCLthreshold); + hist_MSUM_Trigger_Eff = tfs->make(histname4.c_str(),histname4.c_str(), int((MTCAEnd-MTCAStart)/MTCAStep), MTCAStart, MTCAEnd ); + } + if(fMakeTree) + { + //Tree of objects to save + ResetTree(); //call at start of every event + evtTree = tfs->make("TriggerMetrics","Trigger Metric Tree"); + evtTree->Branch("event",&tree_event); + evtTree->Branch("run",&tree_run); + evtTree->Branch("subrun",&tree_subrun); + evtTree->Branch("timestamp",&tree_timestamp); + int MonBins = (fMonStop - fMonStart)/fMonStep+1; + evtTree->Branch("OnBeamMonPeaks",&tree_OnBeamMonPeaks, MonBins, 0); + evtTree->Branch("OnBeamMonPeaksSample",&tree_OnBeamMonPeakSample, MonBins, 0); + evtTree->Branch("OffBeamMonPeaks",&tree_OffBeamMonPeaks, MonBins, 0); + evtTree->Branch("FullBeamMonPeaks",&tree_FullBeamMonPeaks, MonBins, 0); + evtTree->Branch("FullBeamMonPeaksSample",&tree_FullBeamMonPeakSample, MonBins, 0); + evtTree->Branch("MonStart",&tree_MonStart); + evtTree->Branch("MonEnd",&tree_MonEnd); + evtTree->Branch("MonStep",&tree_MonStep); + evtTree->Branch("nAboveThreshold",&tree_nAboveThreshold); + evtTree->Branch("trig_ts",&tree_trig_ts); + evtTree->Branch("promptPE",&tree_promptPE); + evtTree->Branch("prelimPE",&tree_prelimPE); + evtTree->Branch("peakPE",&tree_peakPE); + evtTree->Branch("peaktime",&tree_peaktime); + evtTree->Branch("HLT", &tree_HLTs); + evtTree->Branch("HLTTimes", &tree_HLTTimes); + evtTree->Branch("LLT", &tree_LLTs); + evtTree->Branch("LLTTimes", &tree_LLTTimes); + } + } + + void BeamRateCalib::analyze(art::Event const & e) + { + // Implementation of required member function here. + art::ServiceHandle tfs; //Common art service should read about + ResetTree(); + fEvNumber = e.id().event(); + tree_event = fEvNumber; + tree_run = e.run(); + tree_subrun = e.subRun(); + std::cout << " on run " << tree_run << " event " << fEvNumber << std::endl; + unsigned long long int tsval = e.time().value(); + ETrigTime=0; + const unsigned long int mask32 = 0xFFFFFFFFUL; + tree_timestamp = ( tsval >> 32 ) & mask32; + art::Handle< std::vector< raw::OpDetWaveform > > waveHandle; //User handle for vector of OpDetWaveforms + e.getByLabel(fInputModuleName, fInputInstanceName, fInputProcessName, waveHandle); + if(!waveHandle.isValid() || waveHandle->size() == 0){ + std::cout << "Bad wavehandle found on event " << fEvNumber << " " << waveHandle.isValid() << " " << waveHandle->size() << std::endl; + evtTree->Fill(); + return; + } + art::Handle> ptbHandle; + e.getByLabel(fPTBLabel,ptbHandle); + int HLTIndexer = 0; + for(int index=0; indexsize()); index++) + { + auto ptb = (*ptbHandle)[index]; + auto hltrigs = ptb.GetHLTriggers(); + for(int HLT=0; HLT GrabbedHLTs; + while(Power<32) + { + if(hltrigs[HLT].trigger_word & (0x1 << Power)) + { + GrabbedHLTs.push_back(Power); + if( (Power==1 || Power==2 || Power==3 || Power==4) && (ETrigTime==0 || (hltrigs[HLT].timestamp*20+fEtrigOffset)=int(tree_HLTs.size())) break; //Error handling for weird events + tree_HLTs[HLTIndexer] = GrabbedHLTs[i]; + tree_HLTTimes[HLTIndexer] = hltrigs[HLT].timestamp*20; + HLTIndexer=HLTIndexer+1; + } + } + } + int LLTIndexer=0; + for(int index=0; indexsize()); index++) + { + auto ptb = (*ptbHandle)[index]; + auto lltrigs = ptb.GetLLTriggers(); + for(int LLT=0; LLT GrabbedLLTs; + while(Power<32) + { + if(lltrigs[LLT].trigger_word & (0x1 << Power)) + { + GrabbedLLTs.push_back(Power); + if(Power==26 || Power==30) GateOpenTime=lltrigs[LLT].timestamp*20;//-fCAENOffset; //UNIX ns time of gate opening. Adjusted for PTB-TDC offset + } + Power=Power+1; + } + for(int i=0; i=int(tree_LLTs.size())) continue; //Error handling for weird events + tree_LLTs[LLTIndexer] = GrabbedLLTs[i]; + tree_LLTTimes[LLTIndexer] = lltrigs[LLT].timestamp*20; + LLTIndexer=LLTIndexer+1; + } + } + } + //Seems to be a vector with some extra stuff like isValid() + //raw::OpDetWaveform is a class with start time, vector of samples, and channel ID + analyzeTrigger(waveHandle); //Update the TH2D in other loop...Although maybe we dont need to + if(EventCounter<20 || fSaveAllMON) + { + //SaveChannelWaveforms(waveHandle, fEvNumber); + SaveMonWaveforms(waveHandle, fFCLthreshold, fEvNumber); + } + if(fCheckTriggers) + { + Validate_MonSim_Outputs(e); + } + EventCounter=EventCounter+1; + PMTReadouts = PMTReadouts+(waveHandle->size())/(fTotalCAENBoards*PMTPerBoard); + if(fCheckSoftTrig) + { + art::Handle >softwareTrigHandle; //User handle for vector of OpDetWaveforms + e.getByLabel(fSoftTrigLabel, softwareTrigHandle); + tree_nAboveThreshold = (*softwareTrigHandle)[0].nAboveThreshold; + tree_trig_ts = (*softwareTrigHandle)[0].trig_ts; + tree_promptPE = (*softwareTrigHandle)[0].promptPE; + tree_prelimPE = (*softwareTrigHandle)[0].prelimPE; + tree_peakPE = (*softwareTrigHandle)[0].peakPE; + tree_peaktime = (*softwareTrigHandle)[0].peaktime; + std::cout << "Soft trig size " << (*softwareTrigHandle).size() << " peak PE " << (*softwareTrigHandle)[0].peakPE << " peak time " << (*softwareTrigHandle)[0].peaktime << std::endl; + } + if(fMakeTree) evtTree->Fill(); //call at end of every event + } + + void BeamRateCalib::SaveMonWaveforms(art::Handle< std::vector< raw::OpDetWaveform > > &waveHandle, int MonThreshold, int EventCounter) + { + + // Implement service + art::ServiceHandle tfs; + art::ServiceHandle fTriggerService; + + int TotalFlash = waveHandle->size()/(fTotalCAENBoards*PMTPerBoard); + int GoodFlashIndex=0; + double SmallestTimestamp=0; + for(int FlashCounter=0; FlashCounter *MonPulse = new std::vector(WaveformSize); // add a waveform size getter? + if(!fSaveAllMON) fTriggerService->ConstructMonPulse(*waveHandle, MonThreshold, MonPulse, true, FlashCounter); + else fTriggerService->ConstructMonPulse(*waveHandle, MonThreshold, MonPulse, false, FlashCounter); //Dont save individual channels for now + std::stringstream histname; + if(GoodFlashIndex==FlashCounter) histname << "event_" << EventCounter <<"_Mon"<<"_"<make(histname.str().c_str(), histname.str().c_str(), + MonPulse->size(), 0.0, MonPulse->size()-1); //so this just breaks + for(unsigned int i = 0; i < MonPulse->size(); i++) { + MonHist->SetBinContent(i + 1, (double)(*MonPulse)[i]); //Loop over waveform and set bin content + } + delete MonPulse; + } + } + + void BeamRateCalib::SaveChannelWaveforms(const raw::OpDetWaveform &wvf, int EventCounter, int FlashCounter) + { + std::stringstream histname; + art::ServiceHandle tfs; + int fChNumber = wvf.ChannelNumber(); //Get the channel number for this waveform + histname.str(std::string()); //Resets string stream to nothing + histname << "event_" << EventCounter + << "_opchannel_" << fChNumber + << "_" << opdetType << "_"<make< TH1D >(histname.str().c_str(), histname.str().c_str(), wvf.size(), 0, wvf.size()-1); + for(unsigned int i = 0; i < wvf.size(); i++) { + wvfHist->SetBinContent(i + 1, (double)wvf[i]); //Loop over waveform and set bin content + } + hist_id++; + + } + + void BeamRateCalib::analyzeTrigger(art::Handle< std::vector< raw::OpDetWaveform > > &waveHandle) + { + + // Implement service + art::ServiceHandle tfs; + art::ServiceHandle fTriggerService; + + //Use fill and not set entries for the th2d interface + //Define some MON pulse to do trigger logics on + //Loop over MON thresholds + int TotalFlash = waveHandle->size()/(fTotalCAENBoards*PMTPerBoard); + int GoodFlashIndex=0; + double SmallestTimestamp=0; + bool GoodFlash; + for(int FlashCounter=0; FlashCounter *MonPulse = new std::vector(WaveformSize); + fTriggerService->ConstructMonPulse(*waveHandle, MONThresholds[i], MonPulse, false, FlashCounter); + //Constructed MON pulse now lets compare it to the MTCA requirement + //Flash just before beam gate can offset waveform + //int WaveformOffset = (*waveHandle)[WaveIndex].TimeStamp()*1000 - fNominalGoodStart; //Waveform timestamps are in us // Also delivers time of sample 0 //ns + //Ideally this should be 0 for normal gate opening. For extended waveform we need it to be positive, and our waveform has ealier timestamp + int GateOpenRelativeETrig = -1*int(ETrigTime - GateOpenTime); + int ModifiedWaveStart = (*waveHandle)[WaveIndex].TimeStamp()*1000+fNominalGoodStart; //Adjust timestamp to be actual time of flash trigger + int WaveformOffset = GateOpenRelativeETrig-ModifiedWaveStart; //Measure gate opening time relative to flash trigger + WaveformOffset = WaveformOffset/2; //sample/2 ns + int BeamAcceptanceStart = fBeamWindowStart+WaveformOffset; + int BeamAcceptanceEnd = fBeamWindowEnd+WaveformOffset; + int PeakMon = *std::max_element(MonPulse->begin()+BeamAcceptanceStart, MonPulse->begin()+BeamAcceptanceEnd); + if(MONThresholds[i]==50) + { + std::cout <<"Event " << EventCounter <<" Peak MON " << PeakMon << std::endl; + std::cout << " I check from " << BeamAcceptanceStart << " - " << BeamAcceptanceEnd << std::endl; + } + int PeakMon_NoRange = *std::max_element(MonPulse->begin(), MonPulse->end()); + int PeakOffBeam = *std::max_element(MonPulse->end()-1-(BeamAcceptanceEnd-BeamAcceptanceStart), MonPulse->end()-1); + //Loop over MTCA thresholds + if(GoodFlash) + { + tree_OnBeamMonPeaks[i] = PeakMon; + tree_OnBeamMonPeakSample[i] = std::distance(MonPulse->begin(), std::max_element(MonPulse->begin()+BeamAcceptanceStart, MonPulse->begin()+BeamAcceptanceEnd)); + tree_OffBeamMonPeaks[i] = PeakOffBeam; + tree_FullBeamMonPeakSample[i] = std::distance(MonPulse->begin(), std::max_element(MonPulse->begin(), MonPulse->end()) ); + tree_FullBeamMonPeaks[i] = PeakMon_NoRange; + } + for(int j=0; jsize()); Q++) + { + if( ((*MonPulse)[Q] >=PairVal) && ((*MonPulse)[Q-1] Fill(Q); //Remake tophat plot to look at every crossing + } + } + if(PeakMon>=PairVal && GoodFlash) //Right trigger type here + { + //fill histogram + double FillX = MONThresholds[i]; + double FillY = MTCA_Thresholds[j]; + hist_PropTriggers->Fill(FillX, FillY); //Fill histogram if we trigger + } + if(PeakOffBeam>=PairVal && GoodFlash) //Right trigger type here + { + double FillX = MONThresholds[i]; + double FillY = MTCA_Thresholds[j]; + hist_PropTriggers_OffBeam->Fill(FillX, FillY); + } + if(PeakMon_NoRange>=PairVal && GoodFlash) //Runs over all flashes + { + double FillX = MONThresholds[i]; + double FillY = MTCA_Thresholds[j]; + hist_PropTriggers_NoRestriction->Fill(FillX, FillY); //Fill histogram if we trigger + } + //else // out of things to fill in break for time + //{ + // break; + //} + } //Loop over MTCA threasholds + delete MonPulse; + } //Loop over MON thresholds + } // Loop over flash counters + + } + +//Function made to run over run 15670 where we digitized the MTCA outputs in the timing CAEN +//Assumes we have run the pmttriggerproducer as well +//Validates MON construction in producer, and extracts MTCA threshold settings. Records into histogram + void BeamRateCalib::Validate_MonSim_Outputs(art::Event const & e) + { + /* //Load in waveforms + art::Handle< std::vector< raw::OpDetWaveform > > waveHandle; //User handle for vector of OpDetWaveforms + e.getByLabel(fInputModuleName, fInputInstanceName, fInputProcessName, waveHandle); + //Load in trigger information + art::Handle< std::vector > TriggerHandle; //User handle for vector of OpDetWaveforms + e.getByLabel(std::string("pmttriggerproducer"), std::string(""), std::string("pmtTriggerproducer"), waveHandle); + //Build my MON pulse equivalent + std::vector *MonPulse = new std::vector(WaveformSize); + ConstructMonPulse(waveHandle, MonThreshold, MonPulse, false); + //Load up the relevant waveforms from the timing CAEN + //Have an annoying sample offset and subsampling to handle + //80% of samples come post-trigger and 20% pre-trigger + */ + //sbnd::comm::pmtTrigger has a vector numPassed that says how many PMT passed MON threshold for this (MON) sample + //Also has an a number for the max number of PMT simultaneously passing (basically max of numPassed vector) + //Get channel 903 (I labeled timing CAEN to be 900+ channel ID) + art::Handle< std::vector< raw::OpDetWaveform > > waveHandle; //User handle for vector of OpDetWaveforms + e.getByLabel(fTimingInstanceName, waveHandle); + bool GoodTrigger=false; //could be an vector when we have many flash + bool OneChannel[] = {false, false, false}; + int count=0; + double Energy=0; + int TimingChannelsIncluded=4; + int TimingBoards=1; + int TotalFlash = waveHandle->size()/(TimingBoards*TimingChannelsIncluded); //should match earlier total flash result but different object + for(int FlashCounter=0; FlashCounter907) ) continue; + //Determine if waveform has a pulse in it + std::stringstream histname; + art::ServiceHandle tfs; + int k=0; + //Save waveform + if( (wvf.ChannelNumber()==903 )) SaveChannelWaveforms(wvf, fEvNumber, FlashCounter); + std::vector MonPulse; + while(k < int(wvf.size())) + { + if(wvf[k]-wvf[0] > 100 && TMath::Abs(wvf[wvf.size()-1]-wvf[0])<30) //Make sure no baseline shift + { + if(wvf.ChannelNumber()==905) OneChannel[0] = true; + if(wvf.ChannelNumber()==906) OneChannel[1] = true; + if(wvf.ChannelNumber()==907) OneChannel[2] = true; + //break; + } + if(wvf.ChannelNumber()==903 && wvf[k]-wvf[0] > 17) + { + count=count+1; //We also care about coincident crossings though + } + if(wvf.ChannelNumber()==903) + { + //They disappate more power + //Approximate power through the resistor is ( int(wvf[k]-wvf[0]/30)*0.125 V )**2/50 ohm + //Total energy would be power * 10 us + //Could do both energy and average power + double ADC_Per_Pair = 17; + int wvfMedian = GetWaveformMedian( wvf); + int NPairs = std::max(int(-1*(wvf[k]-wvfMedian)/ADC_Per_Pair), 0); + MonPulse.push_back(NPairs); + double Volts_Per_Pair = 0.125; // 125 mV from CAEN. Knocked down by 10x by attenuators + double Voltage = NPairs*Volts_Per_Pair; + double Resistance = 50; //Ohms + double Power = Voltage*Voltage/Resistance; + Energy =Energy + Power*2*TMath::Power(10, -9); //Power integrated over 2ns. Joules + } + k=k+1; + } + if(wvf.ChannelNumber()==903) + { + double DutyCycle = double(count)/double(wvf.size()); + hist_MSUM_DutyCycle->Fill(DutyCycle); + hist_MSUM_Energy->Fill(Energy); //Can convert to average power by dividing by readout time + //std::cout << Energy << std::endl; + //Make th1d of trigger rate plot of column enabled in fcl + //int WaveformOffset = GateOpenTime - (wvf.TimeStamp()*1000 + ETrigTime)+ fNominalGoodStart ; + //WaveformOffset = WaveformOffset/2; //sample/2 ns + int BeamAcceptanceStart = fBeamWindowStart; // + WaveformOffset; + int BeamAcceptanceEnd = fBeamWindowEnd; // + WaveformOffset; + int PeakMon = *std::max_element(MonPulse.begin()+BeamAcceptanceStart, MonPulse.begin()+BeamAcceptanceEnd); //Grabbing wrong indecies but I am not requiring it to be the trigger pulse. Can add that back later + for(int j=0; j=PairVal) + { + //fill histogram + double FillY = MTCA_Thresholds[j]+float((MTCA_Thresholds[1]-MTCA_Thresholds[0])/2.); + hist_MSUM_Trigger_Eff->Fill(FillY); //Fill histogram if we trigger + } + } + } + }//Channel loop + } //CAEN loop + } //flash loop + + + GoodTrigger = OneChannel[0] && OneChannel[1] && OneChannel[2]; + std::stringstream histname; + art::ServiceHandle tfs; + int FlashCount=0; + if(fTriggerProdCheck) + { + art::Handle< std::vector > TriggerHandle; + e.getByLabel(std::string("pmttriggerproducer"), std::string(""), std::string("PMTDecodeAndTriggerProd"), TriggerHandle); + for(auto& Trig: (*TriggerHandle)) + { + histname.str(std::string()); //Resets string stream to nothing + histname << "event_" << fEvNumber << "pmtriggerprod" << "_flash_" << FlashCount; + std::vector< int > numPassed = Trig.numPassed; + TH1D *pmtTrigHist = tfs->make< TH1D >(histname.str().c_str(), histname.str().c_str(), numPassed.size(), 0, numPassed.size()-1); + for(int i=0; iSetBinContent(i+1, (double)numPassed[i]); + } + FlashCount=FlashCount+1; + } + if(GoodTrigger) + { + for(auto& Trig: (*TriggerHandle)) + { + hist_TriggerThreshold->Fill(Trig.maxPMTs); + } + } + } + //Record how proporition of readout where MSUM is over threshold by a fixed amount for a given run + } + + void BeamRateCalib::endJob() + { + float Scale = 1/float(EventCounter); + //Should do some checks on flash timing to only pick out trigger flash + hist_PropTriggers->Scale(Scale); + hist_PropTriggers->GetXaxis()->SetTitle("MON-Sigma Threshold (adc)"); + hist_PropTriggers->GetYaxis()->SetTitle("MTC/A Threshold (Pairs)"); + hist_PropTriggers->SetEntries(double(EventCounter)); + hist_PropTriggers_OffBeam->Scale(Scale); + hist_PropTriggers_OffBeam->GetXaxis()->SetTitle("MON-Sigma Threshold (adc)"); + hist_PropTriggers_OffBeam->GetYaxis()->SetTitle("MTC/A Threshold (Pairs)"); + hist_PropTriggers_OffBeam->SetEntries(double(EventCounter)); + Scale = 1/float(PMTReadouts); + hist_PropTriggers_NoRestriction->Scale(Scale); + hist_PropTriggers_NoRestriction->GetXaxis()->SetTitle("MON-Sigma Threshold (adc)"); + hist_PropTriggers_NoRestriction->GetYaxis()->SetTitle("MTC/A Threshold (Pairs)"); + hist_PropTriggers_NoRestriction->SetEntries(double(PMTReadouts)); + if(fCheckTriggers) + { + hist_MSUM_Trigger_Eff->Scale(Scale); + hist_MSUM_Trigger_Eff->GetXaxis()->SetTitle("MTC/A Threshold"); + hist_MSUM_Trigger_Eff->GetYaxis()->SetTitle("Proportion Events with HLT"); + hist_MSUM_Trigger_Eff->SetEntries(double(EventCounter)); + } + } + + double BeamRateCalib::GetWaveformMedian(raw::OpDetWaveform wvf) + { + std::sort(wvf.begin(), wvf.end()); + int MedianIndex = int(wvf.size()/2); + return wvf[MedianIndex]; + } + + DEFINE_ART_MODULE(opdet::BeamRateCalib) // Magic line that has to be here + +} diff --git a/sbndcode/OpDetSim/TriggerEmulationService.cc b/sbndcode/OpDetSim/TriggerEmulationService.cc new file mode 100644 index 000000000..eca0f05e5 --- /dev/null +++ b/sbndcode/OpDetSim/TriggerEmulationService.cc @@ -0,0 +1,191 @@ +/** + * @file TriggerEmulationService.cc + * @brief Service implementation associated with TriggerEmulationService.h and TriggerEmulationService_service.cc + * @author Nikki Pallat (palla110@umn.edu) + * @date July 7, 2025 + * @see TriggerEmulationService.h + * @ingroup TriggerEmulationService + * + */ + +#include "TriggerEmulationService.h" + +//namespace lar { +namespace calib { + + TriggerEmulationService::TriggerEmulationService(fhicl::ParameterSet const& pset, art::ActivityRegistry& amp) + : fMonWidth(pset.get("MonWidth", 3)), + fTotalCAENBoards(pset.get("TotalCAENBoards", 10)), + PMTPerBoard(pset.get("PMTPerBoard", 16)), + Baseline(pset.get("Baseline", 14250)), + fMC(pset.get("MC", false)) + {} + + TriggerEmulationService::~TriggerEmulationService() + {} // deconstructor + + void TriggerEmulationService::ConstructMonPulse( + std::vector fWaveforms, + int MonThreshold, + std::vector *MonPulse, + bool Saving, + int FlashCounter, + int *numPairsOverThreshold + ) + { + // Loop over the entries in our waveform vector + // We care about getting the pairing correct + + std::fill(MonPulse->begin(), MonPulse->end(), 0); + + if (fMC) { // monte carlo + if (fWaveforms.empty()) { + std::cout << "Empty waveform vector. Exiting ConstructMonPulse." << std::endl; + return; + } + + std::map channel_to_waveform; + for (const auto& wvf : fWaveforms) + channel_to_waveform[wvf.ChannelNumber()] = &wvf; + + std::vector Pair2 = { 6, 8, 10, 12, 14, 16, 36, 38, 40, 60, 62, 66, 68, 70, 84, 86, 88, 90, 92, 94, 114, 116, 118, 138, 140, 144, 146, 148, 162, 164, 168, 170, 172, 192, 194, 196, 216, 218, 220, 222, 224, 226, 240, 242, 246, 248, 250, 270, 272, 274, 294, 296, 298, 300, 302, 304}; + std::vector Pair1 = { 7, 9, 11, 13, 15, 17, 37, 39, 41, 61, 63, 67, 69, 71, 85, 87, 89, 91, 93, 95, 115, 117, 119, 139, 141, 145, 147, 149, 163, 165, 169, 171, 173, 193, 195, 197, 217, 219, 221, 223, 225, 227, 241, 243, 247, 249, 251, 271, 273, 275, 295, 297, 299, 301, 303, 305}; + std::vector Unpaired = {65, 64, 143, 142, 167, 166, 245, 244}; + std::set used_channels; + + // resize + int ReadoutSize = fWaveforms[0].size(); + MonPulse->assign(ReadoutSize, 0); + + int countPairs = 0; + + for (size_t i = 0; i < Pair1.size(); ++i) { + int ch1 = Pair1[i]; + int ch2 = Pair2[i]; + + // skip if either waveform is missing + if (channel_to_waveform.count(ch1) == 0 || channel_to_waveform.count(ch2) == 0) continue; + // skip if already processed + if (used_channels.count(ch1) || used_channels.count(ch2)) continue; + + const auto& wvf1 = *channel_to_waveform[ch1]; + const auto& wvf2 = *channel_to_waveform[ch2]; + + auto bin1 = ConstructBinaryResponse(wvf1, MonThreshold); + auto bin2 = ConstructBinaryResponse(wvf2, MonThreshold); + + bool pairOverThreshold = false; + + for (int j = 0; j < ReadoutSize; ++j) { + if (bin1[j] || bin2[j]) { + (*MonPulse)[j]++; + pairOverThreshold = true; + } + } + + if (pairOverThreshold) countPairs++; + + used_channels.insert(ch1); + used_channels.insert(ch2); + } + + for (int ch : Unpaired) { // Unpaired channels + if (used_channels.count(ch)) continue; + if (channel_to_waveform.count(ch) == 0) continue; + + const auto& wvf = *channel_to_waveform[ch]; + auto bin = ConstructBinaryResponse(wvf, MonThreshold); + + bool pairOverThreshold = false; + + for (int j = 0; j < ReadoutSize; ++j) { + if (bin[j]) { + (*MonPulse)[j]++; + pairOverThreshold = true; + } + } + + if (pairOverThreshold) countPairs++; + + } + + if (numPairsOverThreshold) *numPairsOverThreshold = countPairs; + + } else { // data + if (fWaveforms.empty()) { + std::cout << "Empty waveform vector. Exiting ConstructMonPulse." << std::endl; + return; + } + + int NumFlash = fWaveforms.size() / (fTotalCAENBoards * PMTPerBoard); + int FirstReadoutIndex = 0 + FlashCounter*PMTPerBoard + 0*PMTPerBoard*NumFlash; + int ReadoutSize = fWaveforms[FirstReadoutIndex].size(); + + for (int CurrentBoard = 0; CurrentBoard < fTotalCAENBoards; ++CurrentBoard) { + int CAENChannel = 0; + // Loop over each PMT in a board + while (CAENChannel < PMTPerBoard) { + + int ChannelStep = 1; + std::vector BinaryMonContrib(ReadoutSize); + + if (CAENChannel != 14) { + int WaveIndex_Pair1 = CAENChannel + FlashCounter*PMTPerBoard + CurrentBoard*PMTPerBoard*NumFlash; + int WaveIndex_Pair2 = CAENChannel + 1 + FlashCounter*PMTPerBoard + CurrentBoard*PMTPerBoard*NumFlash; + + ChannelStep = 2; + + auto const& wvf_Pair1 = fWaveforms[WaveIndex_Pair1]; + auto const& wvf_Pair2 = fWaveforms[WaveIndex_Pair2]; + + std::vector BinaryResponse_Pair1 = ConstructBinaryResponse(wvf_Pair1, MonThreshold); + std::vector BinaryResponse_Pair2 = ConstructBinaryResponse(wvf_Pair2, MonThreshold); + + for (int i=0; i BinaryResponse_Unpaired = ConstructBinaryResponse(wvf_Unpaired, MonThreshold); + for (int i=0; i TriggerEmulationService::ConstructBinaryResponse(const raw::OpDetWaveform &wvf, int MonThreshold) + { + std::vector BinaryResponse(wvf.size()); + int WaveformIndex=1; + while (WaveformIndex-MonThreshold) && ((wvf[WaveformIndex]-Baseline)<-MonThreshold)); + if (CrossedThreshold) + { + int StartIndex = WaveformIndex; + int EndIndex = WaveformIndex+4*fMonWidth; + if (StartIndex>=int(wvf.size() )) StartIndex=wvf.size()-1; //should never happen + if (EndIndex>=int(wvf.size())) EndIndex=wvf.size()-1; //May happen + for (int k=StartIndex; k +#include + +//namespace lar { + namespace calib { + + class TriggerEmulationService { + public: + TriggerEmulationService(fhicl::ParameterSet const&pset, art::ActivityRegistry&); + ~TriggerEmulationService(); + void ConstructMonPulse( + std::vector fWaveforms, + int MonThreshold, + std::vector *MonPulse, + bool Saving, + int FlashCounter, + int *numPairsOverThreshold = nullptr + ); + + int getTotalCAENBoards() const { return fTotalCAENBoards; } + int getPMTPerBoard() const { return PMTPerBoard; } + + //--------------------------------------------------------------------- + private: + std::vector ConstructBinaryResponse(const raw::OpDetWaveform &wvf, int MonThreshold); + + int fMonWidth; + int fTotalCAENBoards; + int PMTPerBoard; + int Baseline; + bool fMC; + }; // TriggerEmulationService + + } // namespace calib + +DECLARE_ART_SERVICE(calib::TriggerEmulationService, LEGACY) +#endif // TRIGGEREMULATION_SERVICE_H + + diff --git a/sbndcode/OpDetSim/TriggerEmulationService_service.cc b/sbndcode/OpDetSim/TriggerEmulationService_service.cc new file mode 100644 index 000000000..28c92d29a --- /dev/null +++ b/sbndcode/OpDetSim/TriggerEmulationService_service.cc @@ -0,0 +1,4 @@ +#include "TriggerEmulationService.h" +#include "art/Framework/Services/Registry/ServiceDefinitionMacros.h" +DEFINE_ART_SERVICE(calib::TriggerEmulationService) + diff --git a/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc b/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc old mode 100644 new mode 100755 index 9e1b6bfd3..f31731954 --- a/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc +++ b/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc @@ -33,6 +33,7 @@ #include #include #include +#include #include "lardataobj/RawData/OpDetWaveform.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" @@ -54,6 +55,8 @@ #include "sbndcode/OpDetSim/opDetSBNDTriggerAlg.hh" #include "sbndcode/OpDetSim/opDetDigitizerWorker.hh" +#include "TriggerEmulationService.h" + namespace opdet { /* @@ -114,6 +117,36 @@ namespace opdet { 1 }; + fhicl::Atom ticksPerSlice{ + Name("ticksPerSlice"), + Comment("Number of ticks for width of sliced waveform. Width of waveform in data is 5000 samples = 10us, since 1 tick of PMT data is 2ns of data.") + }; + + fhicl::Atom PercentTicksBeforeCross{ + Name("PercentTicksBeforeCross"), + Comment("Given a certain width of the waveform, how much of the waveform width is before the crossing point. To match data, this should be 0.2.") + }; + + fhicl::Atom MonThreshold{ + Name("MonThreshold"), + Comment("Threshold") + }; + + fhicl::Atom MonPulseThresh{ + Name("MonPulseThresh"), + Comment("Threshold for MonPulse (to determine interesting trigger)") + }; + + fhicl::Atom SaveNewPlots{ + Name("SaveNewPlots"), + Comment("Save plots of triggered waveforms and MonPulse with new trigger logic.") + }; + + fhicl::Atom SaveOldPlots{ + Name("SaveOldPlots"), + Comment("Save plots of triggered waveforms with old trigger logic.") + }; + fhicl::TableFragment pmtAlgoConfig; fhicl::TableFragment araAlgoConfig; fhicl::TableFragment trigAlgoConfig; @@ -135,13 +168,28 @@ namespace opdet { // Required functions. void produce(art::Event & e) override; - + std::vector sliceWaveforms(std::vector fWaveforms, + int WaveIndex, + std::vector *MonPulse, + int MonPulseThresh, + double tickPeriod, + int ticksPerSlice, + float PercentTicksBeforeCross, + int PMTPerBoard); + std::vector> sliceMonPulse(std::vector *MonPulse, + int MonPulseThresh, + double tickPeriod, + int ticksPerSlice, + float PercentTicksBeforeCross); + void PlotWaveforms(const std::vector& waveforms, + const std::string& basename); opdet::sbndPDMapAlg map; //map for photon detector types unsigned int nChannels = map.size(); std::vector fWaveforms; // holder for un-triggered waveforms private: bool fApplyTriggers; + art::InputTag fInputModuleName; std::unordered_map< raw::Channel_t, std::vector > fFullWaveforms; bool fUseSimPhotonsLite; @@ -152,6 +200,8 @@ namespace opdet { std::vector fWorkers; std::vector> fTriggeredWaveforms; std::vector fWorkerThreads; + std::vector> fSlicedWaveformsAll; + //std::vector> fMonPulsesAll; // product containers std::vector>> fPhotonLiteHandles; @@ -164,6 +214,13 @@ namespace opdet { // trigger algorithm opdet::opDetSBNDTriggerAlg fTriggerAlg; + + int ticksPerSlice; + float PercentTicksBeforeCross; + int MonThreshold; + int MonPulseThresh; + bool SaveNewPlots; + bool SaveOldPlots; }; opDetDigitizerSBND::opDetDigitizerSBND(Parameters const& config) @@ -173,6 +230,12 @@ namespace opdet { , fPMTBaseline(config().pmtAlgoConfig().pmtbaseline()) , fArapucaBaseline(config().araAlgoConfig().baseline()) , fTriggerAlg(config().trigAlgoConfig()) + , ticksPerSlice(config().ticksPerSlice()) + , PercentTicksBeforeCross(config().PercentTicksBeforeCross()) + , MonThreshold(config().MonThreshold()) + , MonPulseThresh(config().MonPulseThresh()) + , SaveNewPlots(config().SaveNewPlots()) + , SaveOldPlots(config().SaveOldPlots()) { opDetDigitizerWorker::Config wConfig( config().pmtAlgoConfig(), config().araAlgoConfig()); @@ -249,6 +312,10 @@ namespace opdet { // Call appropriate produces<>() functions here. produces< std::vector< raw::OpDetWaveform > >(); + produces("triggerEmulation"); + produces("pairsOverThreshold"); + produces< std::vector< raw::OpDetWaveform > >("slicedWaveforms"); + //produces< std::vector< std::vector > >("MonPulses"); } opDetDigitizerSBND::~opDetDigitizerSBND() @@ -294,42 +361,162 @@ namespace opdet { opdet::WaitopDetDigitizerWorkers(fNThreads, fSemFinish); if (fApplyTriggers) { - // find the trigger locations for the waveforms - for (const raw::OpDetWaveform &waveform : fWaveforms) { - raw::Channel_t ch = waveform.ChannelNumber(); - // skip light channels which don't correspond to readout channels - if (ch == std::numeric_limits::max() /* "NULL" value*/) { - continue; - } - raw::ADC_Count_t baseline = (map.isPDType(ch, "pmt_uncoated") || map.isPDType(ch, "pmt_coated")) ? - fPMTBaseline : fArapucaBaseline; - fTriggerAlg.FindTriggerLocations(clockData, detProp, waveform, baseline); - } - // combine the triggers - fTriggerAlg.MergeTriggerLocations(); - // Start the workers! - // Apply the trigger locations - opdet::StartopDetDigitizerWorkers(fNThreads, fSemStart); - opdet::WaitopDetDigitizerWorkers(fNThreads, fSemFinish); - - for (std::vector &waveforms : fTriggeredWaveforms) { - // move these waveforms into the pulseVecPtr - pulseVecPtr->reserve(pulseVecPtr->size() + waveforms.size()); - std::move(waveforms.begin(), waveforms.end(), std::back_inserter(*pulseVecPtr)); - } - // clean up the vector - for (unsigned i = 0; i < fTriggeredWaveforms.size(); i++) { - fTriggeredWaveforms[i] = std::vector(); - } + // clear previous + fSlicedWaveformsAll.clear(); + //fMonPulsesAll.clear(); + // find the trigger locations for the waveforms using the LArService + if (!fWaveforms.empty()) { + // Implement service + art::ServiceHandle tfs; + art::ServiceHandle fTriggerService; + + int PMTPerBoard = fTriggerService->getPMTPerBoard(); + int fTotalCAENBoards = fTriggerService->getTotalCAENBoards(); + + int TotalFlash = fWaveforms.size() / (fTotalCAENBoards * PMTPerBoard); + + int numPairsOverThreshold = 0; + + // Loop through by flash -> compatible with ConstructMonPulse logic + for (int FlashCounter = 0; FlashCounter < TotalFlash; ++FlashCounter) { + int WaveIndex = FlashCounter*PMTPerBoard; + int WaveformSize = fWaveforms[WaveIndex].size(); + + std::vector* MonPulse = new std::vector(WaveformSize, 0); + + int pairThisFlash = 0; + // Send 3ms waveforms to ConstructMonPulse + fTriggerService->ConstructMonPulse(fWaveforms, MonThreshold, MonPulse, false, FlashCounter, &pairThisFlash); + numPairsOverThreshold = numPairsOverThreshold + pairThisFlash; + + double tickPeriod = sampling_rate(clockData); + + std::vector fSlicedWaveforms = sliceWaveforms(fWaveforms, WaveIndex, MonPulse, MonPulseThresh, tickPeriod, ticksPerSlice, PercentTicksBeforeCross, PMTPerBoard); + std::vector> fSlicedMonPulse = sliceMonPulse(MonPulse, MonPulseThresh, tickPeriod, ticksPerSlice, PercentTicksBeforeCross); + + fSlicedWaveformsAll.push_back(std::move(fSlicedWaveforms)); + //fMonPulsesAll.push_back(*MonPulse); + + if (SaveNewPlots) { + // Save histograms + // Sliced waveforms + for (size_t j; j < fSlicedWaveformsAll.size(); ++j) { + std::stringstream plotname2; + plotname2 << "Sliced_waveforms_" << e.id().event() << "_Mon_" << MonThreshold << "_" << FlashCounter << "_slice" << j; + PlotWaveforms(fSlicedWaveformsAll[j], plotname2.str()); + } + // Long MonPulse + std::stringstream histname; + histname << "Long_event_" << e.id().event() << "_Mon_" << MonThreshold << "_" << FlashCounter; + TH1D* MonHist = tfs->make(histname.str().c_str(), histname.str().c_str(), + MonPulse->size(), 0.0, MonPulse->size() - 1); + for (size_t i = 0; i < MonPulse->size(); i++) { + MonHist->SetBinContent(i + 1, (*MonPulse)[i]); + } + // Sliced MonPulse + for (size_t idx = 0; idx < fSlicedMonPulse.size(); ++idx) { + auto const& vec = fSlicedMonPulse[idx]; + std::stringstream histname; + histname << "Sliced_event_" << e.id().event() << "_Mon_" << MonThreshold << "_" << FlashCounter << "_slice" << idx; + + TH1D* MonHist = tfs->make(histname.str().c_str(), histname.str().c_str(), + vec.size(), 0.0, vec.size() - 1); + for (size_t i = 0; i < vec.size(); i++) { + MonHist->SetBinContent(i + 1, vec[i]); + } + } + } + delete MonPulse; + } - // put the waveforms in the event - e.put(std::move(pulseVecPtr)); - // clear out the triggers - fTriggerAlg.ClearTriggerLocations(); + // find the trigger locations for the waveforms - old version, keeping for validation + for (const raw::OpDetWaveform &waveform : fWaveforms) { + raw::Channel_t ch = waveform.ChannelNumber(); + // skip light channels which don't correspond to readout channels + if (ch == std::numeric_limits::max() /* "NULL" value*/) { + continue; + } + raw::ADC_Count_t baseline = (map.isPDType(ch, "pmt_uncoated") || map.isPDType(ch, "pmt_coated")) ? + fPMTBaseline : fArapucaBaseline; + fTriggerAlg.FindTriggerLocations(clockData, detProp, waveform, baseline); + } + + // combine the triggers + fTriggerAlg.MergeTriggerLocations(); + // Start the workers! + // Apply the trigger locations + opdet::StartopDetDigitizerWorkers(fNThreads, fSemStart); + opdet::WaitopDetDigitizerWorkers(fNThreads, fSemFinish); + + // plot fTriggeredWaveforms + if (SaveOldPlots) { + for (size_t j; j < fTriggeredWaveforms.size(); ++j) { + std::stringstream plotnameTW; + plotnameTW << "Triggered_waveforms_" << e.id().event() << "_Mon_" << MonThreshold; + PlotWaveforms(fTriggeredWaveforms[j], plotnameTW.str()); + } + } + + // put triggered waveforms in the event (old trigger logic) + for (std::vector &waveforms : fTriggeredWaveforms) { + // move these waveforms into the pulseVecPtr + pulseVecPtr->reserve(pulseVecPtr->size() + waveforms.size()); + std::move(waveforms.begin(), waveforms.end(), std::back_inserter(*pulseVecPtr)); + } + // clean up the vector + for (unsigned i = 0; i < fTriggeredWaveforms.size(); i++) { + fTriggeredWaveforms[i] = std::vector(); + } + // put the waveforms in the event + e.put(std::move(pulseVecPtr)); + // clear out the triggers + fTriggerAlg.ClearTriggerLocations(); + + + // put boolean trigger result in the event + bool passedTrigger = false; + // passes trigger if any of the fSlicedWaveforms have size > 0 + for (auto wav : fSlicedWaveformsAll) if (wav.size() > 0) passedTrigger = true; + auto triggerFlag = std::make_unique(passedTrigger); + e.put(std::move(triggerFlag), "triggerEmulation"); + + + // put trigger pair result in the event + auto pairFlag = std::make_unique(numPairsOverThreshold); + e.put(std::move(pairFlag), "pairsOverThreshold"); + + + // put sliced waveforms in the event + std::unique_ptr< std::vector< raw::OpDetWaveform > > SlicedWaveformsPtr(std::make_unique< std::vector< raw::OpDetWaveform > > ()); + for (std::vector &waveforms : fSlicedWaveformsAll) { + // move sliced waveforms into the SlicedWaveformsPtr + SlicedWaveformsPtr->reserve(SlicedWaveformsPtr->size() + waveforms.size()); + std::move(waveforms.begin(), waveforms.end(), std::back_inserter(*SlicedWaveformsPtr)); + } + // clean up the vector + for (unsigned i = 0; i < fSlicedWaveformsAll.size(); i++) { + fSlicedWaveformsAll[i] = std::vector(); + } + // put the waveforms in the event + e.put(std::move(SlicedWaveformsPtr), "slicedWaveforms"); + + // put MonPulses in the event + /*auto MonPulsesPtr = std::make_unique>(); + for (auto& pulses : fMonPulsesAll) { + MonPulsesPtr->reserve(MonPulsesPtr->size() + pulses.size()); + std::move(pulses.begin(), pulses.end(), std::back_inserter(*MonPulsesPtr)); + } + // clean up the vector + for (auto& pulses : fMonPulsesAll) pulses.clear(); + // put the waveforms in the event + e.put(std::move(MonPulsesPtr), "MonPulses");*/ + + } else std::cout << "Empty waveforms found on event " << e.id().event() << " " << fWaveforms.empty() << std::endl; } else { + std::cout<<"ApplyTriggers is false"<::max() /* "NULL" value*/) { @@ -345,6 +532,171 @@ namespace opdet { }//produce end + + // sliced MonPulse function: same logic as sliceWaveforms function + std::vector> opDetDigitizerSBND::sliceMonPulse(std::vector *MonPulse, + int MonPulseThresh, + double tickPeriod, + int ticksPerSlice, + float PercentTicksBeforeCross + ) + { + // Slice up each waveform into 10us chunks based on if "interesting" or not + // before and after crossing point (default is ~20% and ~80%) + int ticksBeforeCross = static_cast(std::round(PercentTicksBeforeCross*ticksPerSlice)); + int ticksAfterCross = ticksPerSlice-ticksBeforeCross; + // find interesting area + std::vector> interestIntervals; + std::vector crossingPoints; + // clear + interestIntervals.clear(); + crossingPoints.clear(); + + bool interest = false; + for (int i = 0; i < (int)MonPulse->size(); ++i) { + // find crossing point + if ((*MonPulse)[i] > MonPulseThresh && interest == false) { + crossingPoints.push_back(i); + interest = true; + } else if ((*MonPulse)[i] <= MonPulseThresh) interest = false; + } + + // create 10us slices around crossingPoints + for (int j = 0; j < (int)crossingPoints.size(); ++j) { + + // if near end of full waveform + if (crossingPoints[j]+ticksAfterCross > static_cast(MonPulse->size())) { + interestIntervals.emplace_back(static_cast(MonPulse->size())-ticksPerSlice, static_cast(MonPulse->size())); + } + // if near beginning of full waveform + else if (crossingPoints[j]-ticksBeforeCross < 0) { + interestIntervals.emplace_back(0, ticksPerSlice); + } + else { + // check if overlaps with previous interval + if (!interestIntervals.empty()) { + if (crossingPoints[j]-ticksBeforeCross < interestIntervals.back().second) { + // if overlaps, extend interval + interestIntervals.back() = {interestIntervals.back().first, crossingPoints[j]+ticksAfterCross}; + // if does not overlap, use typical interval length + } else interestIntervals.emplace_back(crossingPoints[j]-ticksBeforeCross, crossingPoints[j]+ticksAfterCross); + // if first, use typical interval length + } else interestIntervals.emplace_back(crossingPoints[j]-ticksBeforeCross, crossingPoints[j]+ticksAfterCross); + } + } + + std::vector> fSlicedMonPulses; + fSlicedMonPulses.clear(); + // loop through intervals + for (auto [start, end] : interestIntervals) { + std::vector sliceMonPulse((*MonPulse).begin() + start, (*MonPulse).begin() + end); + if (!sliceMonPulse.empty()) { + fSlicedMonPulses.push_back(std::move(sliceMonPulse)); + } + } + + return fSlicedMonPulses; + } + + + // sliceWaveforms function + std::vector opDetDigitizerSBND::sliceWaveforms(std::vector fWaveforms, + int WaveIndex, + std::vector *MonPulse, + int MonPulseThresh, + double tickPeriod, + int ticksPerSlice, + float PercentTicksBeforeCross, + int PMTPerBoard + ) + { + // Slice up each waveform into 10us chunks based on if "interesting" or not + // how many ticks correspond to 10us + // before and after crossing point (default is ~20% and ~80%) + int ticksBeforeCross = static_cast(std::round(PercentTicksBeforeCross*ticksPerSlice)); + int ticksAfterCross = ticksPerSlice-ticksBeforeCross; + // find interesting area + std::vector> interestIntervals; + std::vector crossingPoints; + // clear initial variables + crossingPoints.clear(); + interestIntervals.clear(); + + bool interest = false; + for (int i = 0; i < (int)MonPulse->size(); ++i) { + // find crossing point + if ((*MonPulse)[i] > MonPulseThresh && interest == false) { + crossingPoints.push_back(i); + interest = true; + } else if ((*MonPulse)[i] <= MonPulseThresh) interest = false; + } + + // create 10us slices around crossingPoints + for (int j = 0; j < (int)crossingPoints.size(); ++j) { + + // if near end of full waveform + if (crossingPoints[j]+ticksAfterCross > static_cast(MonPulse->size())) { + interestIntervals.emplace_back(static_cast(MonPulse->size())-ticksPerSlice, static_cast(MonPulse->size())); + } + // if near beginning of full waveform + else if (crossingPoints[j]-ticksBeforeCross < 0) { + interestIntervals.emplace_back(0, ticksPerSlice); + } + else { + // check if overlaps with previous interval + if (!interestIntervals.empty()) { + if (crossingPoints[j]-ticksBeforeCross < interestIntervals.back().second) { + // if overlaps, extend interval + interestIntervals.back() = {interestIntervals.back().first, crossingPoints[j]+ticksAfterCross}; + // if does not overlap, use typical interval length + } else interestIntervals.emplace_back(crossingPoints[j]-ticksBeforeCross, crossingPoints[j]+ticksAfterCross); + // if first, use typical interval length + } else interestIntervals.emplace_back(crossingPoints[j]-ticksBeforeCross, crossingPoints[j]+ticksAfterCross); + } + } + + std::vector fSlicedWaveforms; + fSlicedWaveforms.clear(); + // loop through channels + for (int chan = 0; chan < PMTPerBoard; ++chan) { + const raw::OpDetWaveform& wf = fWaveforms[WaveIndex + chan]; + + for (auto [start, end] : interestIntervals) { + double sliceTime = wf.TimeStamp() + start * tickPeriod; + std::vector sliceData(wf.begin() + start, wf.begin() + end); + + if (!sliceData.empty()) { + raw::OpDetWaveform slice(sliceTime, wf.ChannelNumber(), sliceData); + fSlicedWaveforms.push_back(std::move(slice)); + } + } + } + + return fSlicedWaveforms; + } + + void opDetDigitizerSBND::PlotWaveforms(const std::vector& waveforms, + const std::string& basename) + { + art::ServiceHandle tfs; + + for (size_t i = 0; i < waveforms.size(); ++i) { + const auto& wf = waveforms[i]; + + // unique name per event + waveform index + channel + std::stringstream histName; + histName << basename << "_wf" << i << "_chan" << wf.ChannelNumber(); + + TH1F* h = tfs->make(histName.str().c_str(), + Form("OpDet waveform (chan %d, wf %zu)", wf.ChannelNumber(), i), + wf.size(), 0, wf.size()); + + for (size_t tick = 0; tick < wf.size(); ++tick) { + h->SetBinContent(tick + 1, wf[tick]); + } + } + } + DEFINE_ART_MODULE(opdet::opDetDigitizerSBND) }//closing namespace diff --git a/sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl b/sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl old mode 100644 new mode 100755 index 066d83f9e..49da181ad --- a/sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl +++ b/sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl @@ -11,6 +11,13 @@ sbnd_opdetdigitizer: InputModule: "PDFastSim" WaveformSize: 2000.0 #ns (dummy value, resized according to readout of each PD) UseSimPhotonsLite: true # false for SimPhotons + ApplyTriggers: true #optional + ticksPerSlice: 5000 # corresponds to 10us + PercentTicksBeforeCross: 0.2 + MonThreshold: 15 + MonPulseThresh: 15 + SaveNewPlots: true + SaveOldPlots: false @table::sbnd_digipmt_alg @table::sbnd_digiarapuca_alg diff --git a/sbndcode/OpDetSim/run_BeamRatesCalib.fcl b/sbndcode/OpDetSim/run_BeamRatesCalib.fcl new file mode 100644 index 000000000..52cf85b56 --- /dev/null +++ b/sbndcode/OpDetSim/run_BeamRatesCalib.fcl @@ -0,0 +1,37 @@ +#include "detsimmodules.fcl" +#include "BeamRatesCalib_Defaults.fcl" + +process_name: BeamRateCalib + + +# Source is an art-ROOT file containing CAEN waveforms for each beam spill +source: { + module_type: RootInput +} + +services: +{ + TFileService: { fileName: "BeamRateRes.root" } + TriggerEmulationService: { + service_type: TriggerEmulationService + MonWidth: 12 + TotalCAENBoards: 10 + PMTPerBoard: 16 + Baseline: 14250 + MC: false + } +} + + +# Define and schedule filter/producer/analyzers +physics:{ + + analyzers:{ + BeamCalib: @local::BeamRatesCalib + } + reco: [ ] + + ana: [ BeamCalib ] + trigger_paths: [ ] + end_paths: [ ana ] +} diff --git a/sbndcode/OpDetSim/triggeremulationservice_sbnd.fcl b/sbndcode/OpDetSim/triggeremulationservice_sbnd.fcl new file mode 100644 index 000000000..54a431db2 --- /dev/null +++ b/sbndcode/OpDetSim/triggeremulationservice_sbnd.fcl @@ -0,0 +1,15 @@ +BEGIN_PROLOG + +sbnd_triggeremulation_service: { + TFileService: { fileName: "opDetDigitizerRes.root" } + TriggerEmulationService: { + service_type: TriggerEmulationService + MonWidth: 12 + TotalCAENBoards: 10 + PMTPerBoard: 16 + Baseline: 14250 + MC: true + } +} + +END_PROLOG From 8957e5e23614ca8ebfce975bbf504c6a24be869c Mon Sep 17 00:00:00 2001 From: Nikki Pallat Date: Mon, 22 Sep 2025 10:49:49 -0500 Subject: [PATCH 2/7] Update to CMakeLists --- sbndcode/OpDetSim/CMakeLists.txt | 104 ++++++++++++++++++++++++++++++- 1 file changed, 103 insertions(+), 1 deletion(-) mode change 100644 => 100755 sbndcode/OpDetSim/CMakeLists.txt diff --git a/sbndcode/OpDetSim/CMakeLists.txt b/sbndcode/OpDetSim/CMakeLists.txt old mode 100644 new mode 100755 index bb3b51f8a..8d33304a6 --- a/sbndcode/OpDetSim/CMakeLists.txt +++ b/sbndcode/OpDetSim/CMakeLists.txt @@ -26,8 +26,10 @@ art_make_library( SOURCE DigiArapucaSBNDAlg.cc DigiPMTSBNDAlg.cc opDetDigitizerSBND_module.cc + #BeamRateCalibService_service.cc opDetDigitizerWorker.cc opDetSBNDTriggerAlg.cc + TriggerEmulationService.cc LIBRARIES sbndcode_OpDetSim_sbndPDMapAlg_tool larcore::Geometry_Geometry_service @@ -55,9 +57,72 @@ set ( lardata::Utilities lardataobj::RawData lardataobj::RecoBase + lardataobj::RecoBase + lardataobj::headers + lardata::Utilities_LArFFT_service lardataobj::AnalysisBase lardata::DetectorInfoServices_DetectorClocksServiceStandard_service sbndcode_Utilities_SignalShapingServiceSBND_service + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_Common + sbndaq_artdaq_core::sbndaq-artdaq-core_Obj_SBND + # sbndcode_CRT + # sbndcode_CRTUtils + #sbnobj::SBND_CRT + #sbndcode_CRTUtils + sbnobj::SBND_Timing + sbndcode_GeoWrappers + nurandom::RandomUtils_NuRandomService_service + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art::Framework_Services_Optional_RandomNumberGenerator_service + art::Persistency_Common + art::Persistency_Provenance + art::Utilities + canvas::canvas + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + cetlib::cetlib + CLHEP::CLHEP + ROOT::Core + ROOT::Tree +) + +set ( + NOSIGNALSHAPE_LIBRARIES lardata::Utilities_LArFFT_service + larcorealg::Geometry + larcore::Geometry_Geometry_service + lardata::Utilities + lardataalg::DetectorInfo + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art::Persistency_Common + art::Persistency_Provenance + art_root_io::tfile_support + art_root_io::TFileService_service + art::Framework_Services_System_TriggerNamesService_service + art::Utilities + canvas::canvas + messagefacility::MF_MessageLogger + cetlib::cetlib + cetlib_except::cetlib_except + ROOT::Geom + ROOT::Core +) + +set ( + OTHER_LIBRARIES + sbndcode_OpDetSim + larcore::Geometry_Geometry_service + lardataobj::Simulation + lardata::Utilities + lardataobj::RawData + lardataobj::RecoBase + lardataobj::AnalysisBase + lardata::DetectorInfoServices_DetectorClocksServiceStandard_service + #sbndcode_OpDetSim_BeamRateCalibService_service nurandom::RandomUtils_NuRandomService_service art::Framework_Core art::Framework_Principal @@ -68,12 +133,49 @@ set ( cetlib::cetlib CLHEP::CLHEP ROOT::Core + larcorealg::Geometry + larcore::Geometry_Geometry_service + larsim::Simulation + lardataobj::Simulation + larsim::MCCheater_BackTrackerService_service + lardata::Utilities + larevt::Filters + lardataobj::RawData + lardataobj::RecoBase + lardataobj::AnalysisBase + lardataobj::MCBase + larreco::RecoAlg + lardata::RecoObjects + larpandora::LArPandoraInterface + nusimdata::SimulationBase + art::Framework_Services_Registry + art::Persistency_Common + art::Persistency_Provenance + art::Utilities + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + ROOT::Geom + ROOT::XMLIO + ROOT::Gdml + ROOT::Core + sbndcode_RecoUtils ) + + cet_build_plugin(opHitFinderSBND art::module SOURCE opHitFinderSBND_module.cc LIBRARIES ${MODULE_LIBRARIES}) -cet_build_plugin(opDetDigitizerSBND art::module SOURCE opDetDigitizerSBND_module.cc LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(opDetDigitizerSBND art::module SOURCE opDetDigitizerSBND_module.cc LIBRARIES ${MODULE_LIBRARIES} ${OTHER_LIBRARIES}) cet_build_plugin(wvfAna art::module SOURCE wvfAna_module.cc LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(BeamRatesCalib art::module SOURCE BeamRatesCalib_module.cc LIBRARIES ${MODULE_LIBRARIES} ${OTHER_LIBRARIES}) +cet_build_plugin(PMTNoiseCounter art::module SOURCE PMTNoiseCounter_module.cc LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(FlashWithNeutrino art::module SOURCE FlashWithNeutrinoInTPC.cc LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(SummedWaveformProducer art::module SOURCE Zero_Surpressed_PMT_Summer_module.cc LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(TimeStampDumper art::module SOURCE TimeStampDumper_module.cc LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(michelETagger art::module SOURCE michelETagger_module.cc LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(TriggerEmulationService art::service SOURCE TriggerEmulationService_service.cc TriggerEmulationService.cc LIBRARIES ${NOSIGNALSHAPE_LIBRARIES}) +#cet_build_plugin(SimpleBeamRateCalibAnalyzer art::module SOURCE SimpleBeamRateCalibAnalyzer_module.cc LIBRARIES ${OTHER_LIBRARIES}) + install_headers() install_fhicl() From 66c672fda8dbf2f7d4a42b1f886a94c56858e7f3 Mon Sep 17 00:00:00 2001 From: Nikki Pallat Date: Mon, 22 Sep 2025 10:51:28 -0500 Subject: [PATCH 3/7] Another update to CMakeLists --- sbndcode/OpDetSim/CMakeLists.txt | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sbndcode/OpDetSim/CMakeLists.txt b/sbndcode/OpDetSim/CMakeLists.txt index 8d33304a6..8eab90fc1 100755 --- a/sbndcode/OpDetSim/CMakeLists.txt +++ b/sbndcode/OpDetSim/CMakeLists.txt @@ -168,11 +168,11 @@ cet_build_plugin(opHitFinderSBND art::module SOURCE opHitFinderSBND_module.cc LI cet_build_plugin(opDetDigitizerSBND art::module SOURCE opDetDigitizerSBND_module.cc LIBRARIES ${MODULE_LIBRARIES} ${OTHER_LIBRARIES}) cet_build_plugin(wvfAna art::module SOURCE wvfAna_module.cc LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(BeamRatesCalib art::module SOURCE BeamRatesCalib_module.cc LIBRARIES ${MODULE_LIBRARIES} ${OTHER_LIBRARIES}) -cet_build_plugin(PMTNoiseCounter art::module SOURCE PMTNoiseCounter_module.cc LIBRARIES ${MODULE_LIBRARIES}) -cet_build_plugin(FlashWithNeutrino art::module SOURCE FlashWithNeutrinoInTPC.cc LIBRARIES ${MODULE_LIBRARIES}) -cet_build_plugin(SummedWaveformProducer art::module SOURCE Zero_Surpressed_PMT_Summer_module.cc LIBRARIES ${MODULE_LIBRARIES}) -cet_build_plugin(TimeStampDumper art::module SOURCE TimeStampDumper_module.cc LIBRARIES ${MODULE_LIBRARIES}) -cet_build_plugin(michelETagger art::module SOURCE michelETagger_module.cc LIBRARIES ${MODULE_LIBRARIES}) +#cet_build_plugin(PMTNoiseCounter art::module SOURCE PMTNoiseCounter_module.cc LIBRARIES ${MODULE_LIBRARIES}) +#cet_build_plugin(FlashWithNeutrino art::module SOURCE FlashWithNeutrinoInTPC.cc LIBRARIES ${MODULE_LIBRARIES}) +#cet_build_plugin(SummedWaveformProducer art::module SOURCE Zero_Surpressed_PMT_Summer_module.cc LIBRARIES ${MODULE_LIBRARIES}) +#cet_build_plugin(TimeStampDumper art::module SOURCE TimeStampDumper_module.cc LIBRARIES ${MODULE_LIBRARIES}) +#cet_build_plugin(michelETagger art::module SOURCE michelETagger_module.cc LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(TriggerEmulationService art::service SOURCE TriggerEmulationService_service.cc TriggerEmulationService.cc LIBRARIES ${NOSIGNALSHAPE_LIBRARIES}) #cet_build_plugin(SimpleBeamRateCalibAnalyzer art::module SOURCE SimpleBeamRateCalibAnalyzer_module.cc LIBRARIES ${OTHER_LIBRARIES}) From 38015787ff45ab180287b4d06b6e4c6a4ea46f5d Mon Sep 17 00:00:00 2001 From: Nikki Pallat Date: Mon, 22 Sep 2025 17:58:12 -0500 Subject: [PATCH 4/7] Updates to information saved in the event --- .../OpDetSim/opDetDigitizerSBND_module.cc | 27 ++++++++++--------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc b/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc index f31731954..2f355776c 100755 --- a/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc +++ b/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc @@ -201,7 +201,8 @@ namespace opdet { std::vector> fTriggeredWaveforms; std::vector fWorkerThreads; std::vector> fSlicedWaveformsAll; - //std::vector> fMonPulsesAll; + std::vector fMonPulsesFlat; + std::vector pulseSizes; // product containers std::vector>> fPhotonLiteHandles; @@ -315,7 +316,8 @@ namespace opdet { produces("triggerEmulation"); produces("pairsOverThreshold"); produces< std::vector< raw::OpDetWaveform > >("slicedWaveforms"); - //produces< std::vector< std::vector > >("MonPulses"); + produces< std::vector >("MonPulses"); + produces< std::vector >("MonPulseSizes"); } opDetDigitizerSBND::~opDetDigitizerSBND() @@ -364,7 +366,8 @@ namespace opdet { // clear previous fSlicedWaveformsAll.clear(); - //fMonPulsesAll.clear(); + fMonPulsesFlat.clear(); + pulseSizes.clear(); // find the trigger locations for the waveforms using the LArService if (!fWaveforms.empty()) { // Implement service @@ -396,7 +399,8 @@ namespace opdet { std::vector> fSlicedMonPulse = sliceMonPulse(MonPulse, MonPulseThresh, tickPeriod, ticksPerSlice, PercentTicksBeforeCross); fSlicedWaveformsAll.push_back(std::move(fSlicedWaveforms)); - //fMonPulsesAll.push_back(*MonPulse); + fMonPulsesFlat.insert(fMonPulsesFlat.end(), (*MonPulse).begin(), (*MonPulse).end()); + pulseSizes.push_back(MonPulse->size()); if (SaveNewPlots) { // Save histograms @@ -503,15 +507,12 @@ namespace opdet { // put MonPulses in the event - /*auto MonPulsesPtr = std::make_unique>(); - for (auto& pulses : fMonPulsesAll) { - MonPulsesPtr->reserve(MonPulsesPtr->size() + pulses.size()); - std::move(pulses.begin(), pulses.end(), std::back_inserter(*MonPulsesPtr)); - } - // clean up the vector - for (auto& pulses : fMonPulsesAll) pulses.clear(); - // put the waveforms in the event - e.put(std::move(MonPulsesPtr), "MonPulses");*/ + auto flatPtr = std::make_unique>(std::move(fMonPulsesFlat)); + e.put(std::move(flatPtr), "MonPulses"); + + // put pulseSizes in the event + auto sizesPtr = std::make_unique>(std::move(pulseSizes)); + e.put(std::move(sizesPtr), "MonPulseSizes"); } else std::cout << "Empty waveforms found on event " << e.id().event() << " " << fWaveforms.empty() << std::endl; } From b3643702fb26ccf11e19648d99d3a6388dd5345b Mon Sep 17 00:00:00 2001 From: Nikki Pallat Date: Wed, 24 Sep 2025 13:56:42 -0500 Subject: [PATCH 5/7] Corrections based on Jacob's comments on my PR, including default value changes, numPairsOverThreshold calculation, and saving MonPulse in the PMTDecoder --- sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl | 2 +- sbndcode/OpDetSim/BeamRatesCalib_module.cc | 13 +- sbndcode/OpDetSim/TriggerEmulationService.cc | 24 +- sbndcode/OpDetSim/TriggerEmulationService.h | 5 +- .../OpDetSim/opDetDigitizerSBND_module.cc | 280 ++++++++---------- sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl | 2 +- sbndcode/OpDetSim/run_BeamRatesCalib.fcl | 4 +- 7 files changed, 141 insertions(+), 189 deletions(-) diff --git a/sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl b/sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl index 519ad240b..6e6a4bad2 100644 --- a/sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl +++ b/sbndcode/OpDetSim/BeamRatesCalib_Defaults.fcl @@ -33,7 +33,7 @@ BeamRatesCalib: CheckHardwareTriggers: true CheckSoftTrig: true # Number of CAEN boards on for this run - TotalCAENBoards: 10 + TotalCAENBoards: 8 # Assorted labels to get data products out of Art files SoftTrigLabel: "pmtmetricproducer" InputModule: "pmtdecoder" diff --git a/sbndcode/OpDetSim/BeamRatesCalib_module.cc b/sbndcode/OpDetSim/BeamRatesCalib_module.cc index a6ee5bc75..a73fa78b6 100755 --- a/sbndcode/OpDetSim/BeamRatesCalib_module.cc +++ b/sbndcode/OpDetSim/BeamRatesCalib_module.cc @@ -195,8 +195,7 @@ namespace opdet { //OpDet means optical detector fMTCAStart = p.get("MTCAStart", 3); fMTCAStop = p.get("MTCAStop", 64); fMTCAStep = p.get("MTCAStep", 1); - //PMTPerBoard=15; - PMTPerBoard=16; + PMTPerBoard=15; fBeamWindowStart = p.get("BeamWindowStart", 828+680); fBeamWindowEnd = p.get("BeamWindowEnd", 1688+680); fSaveAllMON = p.get("SaveAllMon", false); @@ -205,8 +204,7 @@ namespace opdet { //OpDet means optical detector fOpDetsToPlot = p.get >("OpDetsToPlot"); fCheckSoftTrig = p.get("CheckSoftTrig", false); fSoftTrigLabel = p.get("SoftTrigLabel", "pmtmetricproducer:"); - //fTotalCAENBoards = p.get("TotalCAENBoards", 8); - fTotalCAENBoards = p.get("TotalCAENBoards", 10); + fTotalCAENBoards = p.get("TotalCAENBoards", 8); fPTBLabel = p.get< std::string >("PTBLabel", "ptbdecoder::DECODE"); fTriggerProdCheck = p.get("TriggerProdCheck", false); fNominalWaveformSize = p.get("NominalWaveformSize", 5000); @@ -475,8 +473,9 @@ namespace opdet { //OpDet means optical detector int WaveIndex = FlashCounter*PMTPerBoard; int WaveformSize = (*waveHandle)[WaveIndex].size(); std::vector *MonPulse = new std::vector(WaveformSize); // add a waveform size getter? - if(!fSaveAllMON) fTriggerService->ConstructMonPulse(*waveHandle, MonThreshold, MonPulse, true, FlashCounter); - else fTriggerService->ConstructMonPulse(*waveHandle, MonThreshold, MonPulse, false, FlashCounter); //Dont save individual channels for now + if(!fSaveAllMON) fTriggerService->ConstructMonPulse(*waveHandle, MonThreshold, MonPulse, FlashCounter); + else fTriggerService->ConstructMonPulse(*waveHandle, MonThreshold, MonPulse, FlashCounter); + // Save MonPulses std::stringstream histname; if(GoodFlashIndex==FlashCounter) histname << "event_" << EventCounter <<"_Mon"<<"_"< *MonPulse = new std::vector(WaveformSize); - fTriggerService->ConstructMonPulse(*waveHandle, MONThresholds[i], MonPulse, false, FlashCounter); + fTriggerService->ConstructMonPulse(*waveHandle, MONThresholds[i], MonPulse, FlashCounter); //Constructed MON pulse now lets compare it to the MTCA requirement //Flash just before beam gate can offset waveform //int WaveformOffset = (*waveHandle)[WaveIndex].TimeStamp()*1000 - fNominalGoodStart; //Waveform timestamps are in us // Also delivers time of sample 0 //ns diff --git a/sbndcode/OpDetSim/TriggerEmulationService.cc b/sbndcode/OpDetSim/TriggerEmulationService.cc index eca0f05e5..e09d1b3e6 100644 --- a/sbndcode/OpDetSim/TriggerEmulationService.cc +++ b/sbndcode/OpDetSim/TriggerEmulationService.cc @@ -14,9 +14,9 @@ namespace calib { TriggerEmulationService::TriggerEmulationService(fhicl::ParameterSet const& pset, art::ActivityRegistry& amp) - : fMonWidth(pset.get("MonWidth", 3)), - fTotalCAENBoards(pset.get("TotalCAENBoards", 10)), - PMTPerBoard(pset.get("PMTPerBoard", 16)), + : fMonWidth(pset.get("MonWidth", 12)), + fTotalCAENBoards(pset.get("TotalCAENBoards", 8)), + PMTPerBoard(pset.get("PMTPerBoard", 15)), Baseline(pset.get("Baseline", 14250)), fMC(pset.get("MC", false)) {} @@ -28,7 +28,6 @@ namespace calib { std::vector fWaveforms, int MonThreshold, std::vector *MonPulse, - bool Saving, int FlashCounter, int *numPairsOverThreshold ) @@ -57,8 +56,6 @@ namespace calib { int ReadoutSize = fWaveforms[0].size(); MonPulse->assign(ReadoutSize, 0); - int countPairs = 0; - for (size_t i = 0; i < Pair1.size(); ++i) { int ch1 = Pair1[i]; int ch2 = Pair2[i]; @@ -74,17 +71,12 @@ namespace calib { auto bin1 = ConstructBinaryResponse(wvf1, MonThreshold); auto bin2 = ConstructBinaryResponse(wvf2, MonThreshold); - bool pairOverThreshold = false; - for (int j = 0; j < ReadoutSize; ++j) { if (bin1[j] || bin2[j]) { (*MonPulse)[j]++; - pairOverThreshold = true; } } - if (pairOverThreshold) countPairs++; - used_channels.insert(ch1); used_channels.insert(ch2); } @@ -96,20 +88,14 @@ namespace calib { const auto& wvf = *channel_to_waveform[ch]; auto bin = ConstructBinaryResponse(wvf, MonThreshold); - bool pairOverThreshold = false; - for (int j = 0; j < ReadoutSize; ++j) { if (bin[j]) { (*MonPulse)[j]++; - pairOverThreshold = true; } } - - if (pairOverThreshold) countPairs++; - } - if (numPairsOverThreshold) *numPairsOverThreshold = countPairs; + if (numPairsOverThreshold) *numPairsOverThreshold = *std::max_element(MonPulse->begin(), MonPulse->end()); } else { // data if (fWaveforms.empty()) { @@ -131,7 +117,7 @@ namespace calib { if (CAENChannel != 14) { int WaveIndex_Pair1 = CAENChannel + FlashCounter*PMTPerBoard + CurrentBoard*PMTPerBoard*NumFlash; - int WaveIndex_Pair2 = CAENChannel + 1 + FlashCounter*PMTPerBoard + CurrentBoard*PMTPerBoard*NumFlash; + int WaveIndex_Pair2 = WaveIndex_Pair1 + 1; ChannelStep = 2; diff --git a/sbndcode/OpDetSim/TriggerEmulationService.h b/sbndcode/OpDetSim/TriggerEmulationService.h index 77ea5ab0e..57a925bf9 100644 --- a/sbndcode/OpDetSim/TriggerEmulationService.h +++ b/sbndcode/OpDetSim/TriggerEmulationService.h @@ -12,18 +12,18 @@ #define TRIGGEREMULATION_SERVICE_H -// potential support libraries +// support libraries #include "fhiclcpp/ParameterSet.h" #include "art/Framework/Services/Registry/ServiceMacros.h" #include "art/Framework/Services/Registry/ActivityRegistry.h" #include "art/Framework/Principal/Handle.h" -//#include "canvas/Persistency/Common/Handle.h" #include "lardataobj/RawData/OpDetWaveform.h" #include "larcore/CoreUtils/ServiceUtil.h" // C/C++ standard libraries #include #include +#include //namespace lar { namespace calib { @@ -36,7 +36,6 @@ std::vector fWaveforms, int MonThreshold, std::vector *MonPulse, - bool Saving, int FlashCounter, int *numPairsOverThreshold = nullptr ); diff --git a/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc b/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc index 2f355776c..7ed743f28 100755 --- a/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc +++ b/sbndcode/OpDetSim/opDetDigitizerSBND_module.cc @@ -132,9 +132,9 @@ namespace opdet { Comment("Threshold") }; - fhicl::Atom MonPulseThresh{ - Name("MonPulseThresh"), - Comment("Threshold for MonPulse (to determine interesting trigger)") + fhicl::Atom PairMultiplicityThreshold{ + Name("PairMultiplicityThreshold"), + Comment("Threshold for pair count threshold for event/flash triggers (to determine interesting trigger)") }; fhicl::Atom SaveNewPlots{ @@ -171,14 +171,13 @@ namespace opdet { std::vector sliceWaveforms(std::vector fWaveforms, int WaveIndex, std::vector *MonPulse, - int MonPulseThresh, + int PairMultiplicityThreshold, double tickPeriod, int ticksPerSlice, float PercentTicksBeforeCross, int PMTPerBoard); std::vector> sliceMonPulse(std::vector *MonPulse, - int MonPulseThresh, - double tickPeriod, + int PairMultiplicityThreshold, int ticksPerSlice, float PercentTicksBeforeCross); void PlotWaveforms(const std::vector& waveforms, @@ -200,8 +199,8 @@ namespace opdet { std::vector fWorkers; std::vector> fTriggeredWaveforms; std::vector fWorkerThreads; - std::vector> fSlicedWaveformsAll; - std::vector fMonPulsesFlat; + std::vector> SlicedWaveformsAll; + std::vector MonPulsesFlat; std::vector pulseSizes; // product containers @@ -219,7 +218,7 @@ namespace opdet { int ticksPerSlice; float PercentTicksBeforeCross; int MonThreshold; - int MonPulseThresh; + int PairMultiplicityThreshold; bool SaveNewPlots; bool SaveOldPlots; }; @@ -234,7 +233,7 @@ namespace opdet { , ticksPerSlice(config().ticksPerSlice()) , PercentTicksBeforeCross(config().PercentTicksBeforeCross()) , MonThreshold(config().MonThreshold()) - , MonPulseThresh(config().MonPulseThresh()) + , PairMultiplicityThreshold(config().PairMultiplicityThreshold()) , SaveNewPlots(config().SaveNewPlots()) , SaveOldPlots(config().SaveOldPlots()) { @@ -365,8 +364,8 @@ namespace opdet { if (fApplyTriggers) { // clear previous - fSlicedWaveformsAll.clear(); - fMonPulsesFlat.clear(); + SlicedWaveformsAll.clear(); + MonPulsesFlat.clear(); pulseSizes.clear(); // find the trigger locations for the waveforms using the LArService if (!fWaveforms.empty()) { @@ -390,25 +389,25 @@ namespace opdet { int pairThisFlash = 0; // Send 3ms waveforms to ConstructMonPulse - fTriggerService->ConstructMonPulse(fWaveforms, MonThreshold, MonPulse, false, FlashCounter, &pairThisFlash); + fTriggerService->ConstructMonPulse(fWaveforms, MonThreshold, MonPulse, FlashCounter, &pairThisFlash); numPairsOverThreshold = numPairsOverThreshold + pairThisFlash; double tickPeriod = sampling_rate(clockData); - std::vector fSlicedWaveforms = sliceWaveforms(fWaveforms, WaveIndex, MonPulse, MonPulseThresh, tickPeriod, ticksPerSlice, PercentTicksBeforeCross, PMTPerBoard); - std::vector> fSlicedMonPulse = sliceMonPulse(MonPulse, MonPulseThresh, tickPeriod, ticksPerSlice, PercentTicksBeforeCross); + std::vector SlicedWaveforms = sliceWaveforms(fWaveforms, WaveIndex, MonPulse, PairMultiplicityThreshold, tickPeriod, ticksPerSlice, PercentTicksBeforeCross, PMTPerBoard); + std::vector> SlicedMonPulse = sliceMonPulse(MonPulse, PairMultiplicityThreshold, ticksPerSlice, PercentTicksBeforeCross); - fSlicedWaveformsAll.push_back(std::move(fSlicedWaveforms)); - fMonPulsesFlat.insert(fMonPulsesFlat.end(), (*MonPulse).begin(), (*MonPulse).end()); + SlicedWaveformsAll.push_back(std::move(SlicedWaveforms)); + MonPulsesFlat.insert(MonPulsesFlat.end(), (*MonPulse).begin(), (*MonPulse).end()); pulseSizes.push_back(MonPulse->size()); if (SaveNewPlots) { // Save histograms // Sliced waveforms - for (size_t j; j < fSlicedWaveformsAll.size(); ++j) { + for (size_t j; j < SlicedWaveformsAll.size(); ++j) { std::stringstream plotname2; plotname2 << "Sliced_waveforms_" << e.id().event() << "_Mon_" << MonThreshold << "_" << FlashCounter << "_slice" << j; - PlotWaveforms(fSlicedWaveformsAll[j], plotname2.str()); + PlotWaveforms(SlicedWaveformsAll[j], plotname2.str()); } // Long MonPulse std::stringstream histname; @@ -419,8 +418,8 @@ namespace opdet { MonHist->SetBinContent(i + 1, (*MonPulse)[i]); } // Sliced MonPulse - for (size_t idx = 0; idx < fSlicedMonPulse.size(); ++idx) { - auto const& vec = fSlicedMonPulse[idx]; + for (size_t idx = 0; idx < SlicedMonPulse.size(); ++idx) { + auto const& vec = SlicedMonPulse[idx]; std::stringstream histname; histname << "Sliced_event_" << e.id().event() << "_Mon_" << MonThreshold << "_" << FlashCounter << "_slice" << idx; @@ -480,8 +479,8 @@ namespace opdet { // put boolean trigger result in the event bool passedTrigger = false; - // passes trigger if any of the fSlicedWaveforms have size > 0 - for (auto wav : fSlicedWaveformsAll) if (wav.size() > 0) passedTrigger = true; + // passes trigger if any of the SlicedWaveforms have size > 0 + for (auto wav : SlicedWaveformsAll) if (wav.size() > 0) passedTrigger = true; auto triggerFlag = std::make_unique(passedTrigger); e.put(std::move(triggerFlag), "triggerEmulation"); @@ -493,21 +492,21 @@ namespace opdet { // put sliced waveforms in the event std::unique_ptr< std::vector< raw::OpDetWaveform > > SlicedWaveformsPtr(std::make_unique< std::vector< raw::OpDetWaveform > > ()); - for (std::vector &waveforms : fSlicedWaveformsAll) { + for (std::vector &waveforms : SlicedWaveformsAll) { // move sliced waveforms into the SlicedWaveformsPtr SlicedWaveformsPtr->reserve(SlicedWaveformsPtr->size() + waveforms.size()); std::move(waveforms.begin(), waveforms.end(), std::back_inserter(*SlicedWaveformsPtr)); } // clean up the vector - for (unsigned i = 0; i < fSlicedWaveformsAll.size(); i++) { - fSlicedWaveformsAll[i] = std::vector(); + for (unsigned i = 0; i < SlicedWaveformsAll.size(); i++) { + SlicedWaveformsAll[i] = std::vector(); } // put the waveforms in the event e.put(std::move(SlicedWaveformsPtr), "slicedWaveforms"); // put MonPulses in the event - auto flatPtr = std::make_unique>(std::move(fMonPulsesFlat)); + auto flatPtr = std::make_unique>(std::move(MonPulsesFlat)); e.put(std::move(flatPtr), "MonPulses"); // put pulseSizes in the event @@ -534,146 +533,115 @@ namespace opdet { }//produce end - // sliced MonPulse function: same logic as sliceWaveforms function - std::vector> opDetDigitizerSBND::sliceMonPulse(std::vector *MonPulse, - int MonPulseThresh, - double tickPeriod, - int ticksPerSlice, - float PercentTicksBeforeCross - ) - { - // Slice up each waveform into 10us chunks based on if "interesting" or not - // before and after crossing point (default is ~20% and ~80%) - int ticksBeforeCross = static_cast(std::round(PercentTicksBeforeCross*ticksPerSlice)); - int ticksAfterCross = ticksPerSlice-ticksBeforeCross; - // find interesting area - std::vector> interestIntervals; - std::vector crossingPoints; - // clear - interestIntervals.clear(); - crossingPoints.clear(); - - bool interest = false; - for (int i = 0; i < (int)MonPulse->size(); ++i) { - // find crossing point - if ((*MonPulse)[i] > MonPulseThresh && interest == false) { - crossingPoints.push_back(i); - interest = true; - } else if ((*MonPulse)[i] <= MonPulseThresh) interest = false; - } - - // create 10us slices around crossingPoints - for (int j = 0; j < (int)crossingPoints.size(); ++j) { - - // if near end of full waveform - if (crossingPoints[j]+ticksAfterCross > static_cast(MonPulse->size())) { - interestIntervals.emplace_back(static_cast(MonPulse->size())-ticksPerSlice, static_cast(MonPulse->size())); - } - // if near beginning of full waveform - else if (crossingPoints[j]-ticksBeforeCross < 0) { - interestIntervals.emplace_back(0, ticksPerSlice); - } - else { - // check if overlaps with previous interval - if (!interestIntervals.empty()) { - if (crossingPoints[j]-ticksBeforeCross < interestIntervals.back().second) { - // if overlaps, extend interval - interestIntervals.back() = {interestIntervals.back().first, crossingPoints[j]+ticksAfterCross}; - // if does not overlap, use typical interval length - } else interestIntervals.emplace_back(crossingPoints[j]-ticksBeforeCross, crossingPoints[j]+ticksAfterCross); - // if first, use typical interval length - } else interestIntervals.emplace_back(crossingPoints[j]-ticksBeforeCross, crossingPoints[j]+ticksAfterCross); - } - } + template + std::vector> + findInterestIntervals(const PulseContainer* pulse, + int PairMultiplicityThreshold, + int ticksBeforeCross, + int ticksAfterCross) + { + // find interesting area + std::vector> interestIntervals; + std::vector crossingPoints; + // clear initial variables + crossingPoints.clear(); + interestIntervals.clear(); + bool interest = false; + + for (int i = 0; i < static_cast(pulse->size()); ++i) { + // find crossing point + if ((*pulse)[i] > PairMultiplicityThreshold && !interest) { + crossingPoints.push_back(i); + interest = true; + } + else if ((*pulse)[i] <= PairMultiplicityThreshold) { + interest = false; + } + } - std::vector> fSlicedMonPulses; - fSlicedMonPulses.clear(); - // loop through intervals - for (auto [start, end] : interestIntervals) { - std::vector sliceMonPulse((*MonPulse).begin() + start, (*MonPulse).begin() + end); - if (!sliceMonPulse.empty()) { - fSlicedMonPulses.push_back(std::move(sliceMonPulse)); - } + // create 10us (or given) slices around crossingPoints + for (int j = 0; j < static_cast(crossingPoints.size()); ++j) { + // if near end of full waveform + if (crossingPoints[j] + ticksAfterCross > static_cast(pulse->size())) { + interestIntervals.emplace_back(crossingPoints[j] - ticksBeforeCross, static_cast(pulse->size()) - 1); + } + // if near beginning of full waveform + else if (crossingPoints[j] - ticksBeforeCross < 0) { + interestIntervals.emplace_back(0, crossingPoints[j] + ticksAfterCross); + } + else { + // check if overlaps with previous interval + if (!interestIntervals.empty() && crossingPoints[j] - ticksBeforeCross < interestIntervals.back().second) { + // if overlaps, extend interval + interestIntervals.back() = {interestIntervals.back().first, crossingPoints[j] + ticksAfterCross}; + } + // if does not overlap or if first, use typical interval length + else { + interestIntervals.emplace_back(crossingPoints[j] - ticksBeforeCross, crossingPoints[j] + ticksAfterCross); } + } + } - return fSlicedMonPulses; + return interestIntervals; } + // sliced MonPulses + std::vector> opDetDigitizerSBND::sliceMonPulse( + std::vector* MonPulse, + int PairMultiplicityThreshold, + int ticksPerSlice, + float PercentTicksBeforeCross) + { + // before and after crossing point (default is ~20% and ~80%) + int ticksBeforeCross = static_cast(std::round(PercentTicksBeforeCross*ticksPerSlice)); + int ticksAfterCross = ticksPerSlice - ticksBeforeCross; + + // Slice up each waveform into 10us (or given) chunks based on if "interesting" or not + auto intervals = findInterestIntervals(MonPulse, PairMultiplicityThreshold, ticksBeforeCross, ticksAfterCross); + + std::vector> SlicedMonPulses; + for (auto [start, end] : intervals) { + std::vector slicedMonPulse(MonPulse->begin() + start, MonPulse->begin() + end); + if (!slicedMonPulse.empty()) SlicedMonPulses.push_back(std::move(slicedMonPulse)); + } + return SlicedMonPulses; + } + // sliceWaveforms function - std::vector opDetDigitizerSBND::sliceWaveforms(std::vector fWaveforms, - int WaveIndex, - std::vector *MonPulse, - int MonPulseThresh, - double tickPeriod, - int ticksPerSlice, - float PercentTicksBeforeCross, - int PMTPerBoard - ) - { - // Slice up each waveform into 10us chunks based on if "interesting" or not - // how many ticks correspond to 10us - // before and after crossing point (default is ~20% and ~80%) - int ticksBeforeCross = static_cast(std::round(PercentTicksBeforeCross*ticksPerSlice)); - int ticksAfterCross = ticksPerSlice-ticksBeforeCross; - // find interesting area - std::vector> interestIntervals; - std::vector crossingPoints; - // clear initial variables - crossingPoints.clear(); - interestIntervals.clear(); - - bool interest = false; - for (int i = 0; i < (int)MonPulse->size(); ++i) { - // find crossing point - if ((*MonPulse)[i] > MonPulseThresh && interest == false) { - crossingPoints.push_back(i); - interest = true; - } else if ((*MonPulse)[i] <= MonPulseThresh) interest = false; - } - - // create 10us slices around crossingPoints - for (int j = 0; j < (int)crossingPoints.size(); ++j) { - - // if near end of full waveform - if (crossingPoints[j]+ticksAfterCross > static_cast(MonPulse->size())) { - interestIntervals.emplace_back(static_cast(MonPulse->size())-ticksPerSlice, static_cast(MonPulse->size())); - } - // if near beginning of full waveform - else if (crossingPoints[j]-ticksBeforeCross < 0) { - interestIntervals.emplace_back(0, ticksPerSlice); - } - else { - // check if overlaps with previous interval - if (!interestIntervals.empty()) { - if (crossingPoints[j]-ticksBeforeCross < interestIntervals.back().second) { - // if overlaps, extend interval - interestIntervals.back() = {interestIntervals.back().first, crossingPoints[j]+ticksAfterCross}; - // if does not overlap, use typical interval length - } else interestIntervals.emplace_back(crossingPoints[j]-ticksBeforeCross, crossingPoints[j]+ticksAfterCross); - // if first, use typical interval length - } else interestIntervals.emplace_back(crossingPoints[j]-ticksBeforeCross, crossingPoints[j]+ticksAfterCross); - } - } + std::vector opDetDigitizerSBND::sliceWaveforms( + std::vector fWaveforms, + int WaveIndex, + std::vector* MonPulse, + int PairMultiplicityThreshold, + double tickPeriod, + int ticksPerSlice, + float PercentTicksBeforeCross, + int PMTPerBoard) + { + // before and after crossing point (default is ~20% and ~80%) + int ticksBeforeCross = static_cast(std::round(PercentTicksBeforeCross*ticksPerSlice)); + int ticksAfterCross = ticksPerSlice - ticksBeforeCross; - std::vector fSlicedWaveforms; - fSlicedWaveforms.clear(); - // loop through channels - for (int chan = 0; chan < PMTPerBoard; ++chan) { - const raw::OpDetWaveform& wf = fWaveforms[WaveIndex + chan]; + // Slice up each waveform into 10us (or given) chunks based on if "interesting" or not + auto intervals = findInterestIntervals(MonPulse, PairMultiplicityThreshold, ticksBeforeCross, ticksAfterCross); - for (auto [start, end] : interestIntervals) { - double sliceTime = wf.TimeStamp() + start * tickPeriod; - std::vector sliceData(wf.begin() + start, wf.begin() + end); + std::vector SlicedWaveforms; + // loop through channels + for (int chan = 0; chan < PMTPerBoard; ++chan) { + const raw::OpDetWaveform& wf = fWaveforms[WaveIndex + chan]; - if (!sliceData.empty()) { - raw::OpDetWaveform slice(sliceTime, wf.ChannelNumber(), sliceData); - fSlicedWaveforms.push_back(std::move(slice)); - } - } + for (auto [start, end] : intervals) { + double sliceTime = wf.TimeStamp() + start * tickPeriod; + std::vector sliceData(wf.begin() + start, wf.begin() + end); + + if (!sliceData.empty()) { + raw::OpDetWaveform slice(sliceTime, wf.ChannelNumber(), sliceData); + SlicedWaveforms.push_back(std::move(slice)); } - - return fSlicedWaveforms; + } + } + return SlicedWaveforms; } void opDetDigitizerSBND::PlotWaveforms(const std::vector& waveforms, diff --git a/sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl b/sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl index 49da181ad..793231f5d 100755 --- a/sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl +++ b/sbndcode/OpDetSim/opdetdigitizer_sbnd.fcl @@ -15,7 +15,7 @@ sbnd_opdetdigitizer: ticksPerSlice: 5000 # corresponds to 10us PercentTicksBeforeCross: 0.2 MonThreshold: 15 - MonPulseThresh: 15 + PairMultiplicityThreshold: 15 SaveNewPlots: true SaveOldPlots: false diff --git a/sbndcode/OpDetSim/run_BeamRatesCalib.fcl b/sbndcode/OpDetSim/run_BeamRatesCalib.fcl index 52cf85b56..39ce40d1a 100644 --- a/sbndcode/OpDetSim/run_BeamRatesCalib.fcl +++ b/sbndcode/OpDetSim/run_BeamRatesCalib.fcl @@ -15,8 +15,8 @@ services: TriggerEmulationService: { service_type: TriggerEmulationService MonWidth: 12 - TotalCAENBoards: 10 - PMTPerBoard: 16 + TotalCAENBoards: 8 + PMTPerBoard: 15 Baseline: 14250 MC: false } From 63d191d6bd72384ca28a8d7ed1d6a44792755d80 Mon Sep 17 00:00:00 2001 From: Nikki Pallat Date: Wed, 24 Sep 2025 14:52:22 -0500 Subject: [PATCH 6/7] verified the decoders can still run over raw data with the new changes --- sbndcode/Decoders/PMT/CMakeLists.txt | 1 + .../Decoders/PMT/SBNDPMTDecoder_module.cc | 54 +++++++++++++++++++ sbndcode/Decoders/PMT/pmtdecoder.fcl | 5 +- sbndcode/Decoders/PMT/run_pmtdecoder.fcl | 10 +++- 4 files changed, 68 insertions(+), 2 deletions(-) diff --git a/sbndcode/Decoders/PMT/CMakeLists.txt b/sbndcode/Decoders/PMT/CMakeLists.txt index 5897e0e87..02ea063b9 100644 --- a/sbndcode/Decoders/PMT/CMakeLists.txt +++ b/sbndcode/Decoders/PMT/CMakeLists.txt @@ -3,6 +3,7 @@ cet_build_plugin( SBNDPMTDecoder art::module LIBRARIES sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_SBND sbnobj::SBND_Timing + sbndcode_OpDetSim lardata::Utilities canvas::canvas lardataobj::RecoBase diff --git a/sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc b/sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc index eeb831509..6f82c8347 100644 --- a/sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc +++ b/sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc @@ -49,6 +49,8 @@ #include #include +#include "sbndcode/OpDetSim/TriggerEmulationService.h" + namespace sbndaq { class SBNDPMTDecoder; } @@ -123,6 +125,8 @@ class sbndaq::SBNDPMTDecoder : public art::EDProducer { std::vector fch_map; + int fmon_threshold; + // histogram info std::stringstream histname; //raw waveform hist name art::ServiceHandle tfs; @@ -177,6 +181,8 @@ sbndaq::SBNDPMTDecoder::SBNDPMTDecoder(fhicl::ParameterSet const& p) fch_map = p.get>("ch_map",{}); + fmon_threshold = p.get("mon_threshold", 15); + produces< std::vector< raw::OpDetWaveform > >(fpmt_instance_name); produces< std::vector< raw::OpDetWaveform > >(fflt_instance_name); produces< std::vector< raw::OpDetWaveform > >(ftim_instance_name); @@ -186,6 +192,9 @@ sbndaq::SBNDPMTDecoder::SBNDPMTDecoder(fhicl::ParameterSet const& p) produces< art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform > >(fpmt_timing_instance_name); produces< art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform > >(fflt_timing_instance_name); produces< art::Assns< raw::pmt::BoardTimingInfo, raw::OpDetWaveform > >(ftim_timing_instance_name); + + produces< std::vector >("MonPulses"); + produces< std::vector >("MonPulseSizes"); } void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) @@ -264,6 +273,11 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) evt.put(std::move(pmtTimingAssns),fpmt_timing_instance_name); evt.put(std::move(fltTimingAssns),fflt_timing_instance_name); evt.put(std::move(timTimingAssns),ftim_timing_instance_name); + + auto flatPtr = std::make_unique>(); + auto sizesPtr = std::make_unique>(); + evt.put(std::move(flatPtr), "MonPulses"); + evt.put(std::move(sizesPtr), "MonPulseSizes"); return; } @@ -597,6 +611,43 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) } } } // end board loop + + // loop through flashes + art::ServiceHandle tfs; + art::ServiceHandle fTriggerService; + int PMTPerBoard = fTriggerService->getPMTPerBoard(); + // int fTotalCAENBoards = fTriggerService->getTotalCAENBoards(); + //std::vector< std::vector > MonPulsesAll; + //MonPulseAll.clear(); + std::vector MonPulsesFlat; + std::vector pulseSizes; + MonPulsesFlat.clear(); + pulseSizes.clear(); + int TotalFlash = pmtwvfmVec->size()/((int)fn_caenboards*PMTPerBoard); // pmtwvfmVec = waveHandle ??? + for (int FlashCounter=0; FlashCounter *MonPulse = new std::vector(WaveformSize); + fTriggerService->ConstructMonPulse(*pmtwvfmVec, fmon_threshold, MonPulse, FlashCounter); + //MonPulsesAll.push_back(std::move(MonPulse)); + MonPulsesFlat.insert(MonPulsesFlat.end(), (*MonPulse).begin(), (*MonPulse).end()); + pulseSizes.push_back(MonPulse->size()); + delete MonPulse; + } + // make ptrs + auto flatPtr = std::make_unique>(std::move(MonPulsesFlat)); + auto sizesPtr = std::make_unique>(std::move(pulseSizes)); + + /*std::unique_ptr< std::vector< std::vector > > MonPulsesPtr(std::make_unique< std::vector< std::vector > > ()); + for (auto &pulse : MonPulsesAll) { + MonPulsesPtr->reserve(MonPulsesPtr->size() + pulse.size()); + std::move(pulse.begin(), pulse.end(), std::back_inserter(*MonPulsesPtr)); + } + // clean up the vector + for (unsigned i = 0; i < MonPulsesAll.size(); i++) MonPulsesAll[i] = std::vector(); + */ + board_frag_v.clear(); evt.put(std::move(pmtwvfmVec),fpmt_instance_name); @@ -609,6 +660,9 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) evt.put(std::move(pmtTimingAssns),fpmt_timing_instance_name); evt.put(std::move(fltTimingAssns),fflt_timing_instance_name); evt.put(std::move(timTimingAssns),ftim_timing_instance_name); + + evt.put(std::move(flatPtr), "MonPulses"); + evt.put(std::move(sizesPtr), "MonPulseSizes"); } void sbndaq::SBNDPMTDecoder::get_fragments(artdaq::Fragment & frag, std::vector> & board_frag_v){ diff --git a/sbndcode/Decoders/PMT/pmtdecoder.fcl b/sbndcode/Decoders/PMT/pmtdecoder.fcl index fbc589a9e..421fe92a5 100644 --- a/sbndcode/Decoders/PMT/pmtdecoder.fcl +++ b/sbndcode/Decoders/PMT/pmtdecoder.fcl @@ -47,6 +47,9 @@ pmtdecoder: fragid_offset: 40960 # the offset to subtract to get single digit fragids; if use_set_map is true, should set this to 0 hist_evt: 1 # the # of the event used to generate the histograms, 1st event by default + # trigger configurable + mon_threshold: 15 # ADC channel value threshold to add 1 to the trigger response MON pulse + # for when the fragIDs are mapped to the old configuration or you need to hardcode fragids... ## to use this, must set `fragid_offset` to 0!!! use_set_map: false @@ -66,4 +69,4 @@ pmtdecoder: 900,901,902,903,904,905,906,907,908,909,910,911,912,913,914] # digitizer 8 (TIMING CAEN) } -END_PROLOG \ No newline at end of file +END_PROLOG diff --git a/sbndcode/Decoders/PMT/run_pmtdecoder.fcl b/sbndcode/Decoders/PMT/run_pmtdecoder.fcl index 844551ab5..aa67d4ae4 100644 --- a/sbndcode/Decoders/PMT/run_pmtdecoder.fcl +++ b/sbndcode/Decoders/PMT/run_pmtdecoder.fcl @@ -6,6 +6,14 @@ process_name: PMTDecoder services: { TFileService : {fileName: "decoder_hist.root"} + TriggerEmulationService: { + service_type: TriggerEmulationService + MonWidth: 12 + TotalCAENBoards: 8 + PMTPerBoard: 15 + Baseline: 14250 + MC: false + } } #source is a root file @@ -51,4 +59,4 @@ physics: # 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 +} From 584880a9b49a6cb9be1c4d7bf64852cad42b7c2d Mon Sep 17 00:00:00 2001 From: Nikki Pallat Date: Fri, 3 Oct 2025 14:50:51 -0500 Subject: [PATCH 7/7] Bug fix: for incorrectly nested services fcl and updates to default fcl parameters --- sbndcode/OpDetSim/TriggerEmulationService.cc | 2 +- .../OpDetSim/triggeremulationservice_sbnd.fcl | 16 +++++++--------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/sbndcode/OpDetSim/TriggerEmulationService.cc b/sbndcode/OpDetSim/TriggerEmulationService.cc index e09d1b3e6..95e24206b 100644 --- a/sbndcode/OpDetSim/TriggerEmulationService.cc +++ b/sbndcode/OpDetSim/TriggerEmulationService.cc @@ -18,7 +18,7 @@ namespace calib { fTotalCAENBoards(pset.get("TotalCAENBoards", 8)), PMTPerBoard(pset.get("PMTPerBoard", 15)), Baseline(pset.get("Baseline", 14250)), - fMC(pset.get("MC", false)) + fMC(pset.get("MC", true)) {} TriggerEmulationService::~TriggerEmulationService() diff --git a/sbndcode/OpDetSim/triggeremulationservice_sbnd.fcl b/sbndcode/OpDetSim/triggeremulationservice_sbnd.fcl index 54a431db2..20c1a3428 100644 --- a/sbndcode/OpDetSim/triggeremulationservice_sbnd.fcl +++ b/sbndcode/OpDetSim/triggeremulationservice_sbnd.fcl @@ -1,15 +1,13 @@ BEGIN_PROLOG sbnd_triggeremulation_service: { - TFileService: { fileName: "opDetDigitizerRes.root" } - TriggerEmulationService: { - service_type: TriggerEmulationService - MonWidth: 12 - TotalCAENBoards: 10 - PMTPerBoard: 16 - Baseline: 14250 - MC: true - } + service_type: TriggerEmulationService + TFileService: { fileName: "opDetDigitizerRes.root" } + MonWidth: 12 + TotalCAENBoards: 8 + PMTPerBoard: 15 + Baseline: 14250 + MC: true } END_PROLOG