From eb71091a7f21019038959a00a4e6abc5596b472d Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Thu, 21 Nov 2024 13:12:56 -0600 Subject: [PATCH 1/9] First commit and it's a big one --- DataModel/Waveform.h | 54 +++ UserTools/Factory/Factory.cpp | 174 +-------- UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 343 ++++++++++++++++++ UserTools/PMTWaveformSim/PMTWaveformSim.h | 75 ++++ UserTools/PMTWaveformSim/README.md | 20 + .../PhaseIIADCCalibrator.cpp | 34 +- .../PhaseIIADCHitFinder.cpp | 82 ++++- .../PhaseIIADCHitFinder/PhaseIIADCHitFinder.h | 3 +- UserTools/Unity.h | 179 +-------- configfiles/PMTWaveformSim/DummyToolConfig | 3 + configfiles/PMTWaveformSim/LoadGeometryConfig | 9 + configfiles/PMTWaveformSim/LoadWCSimConfig | 16 + .../PMTWaveformSim/PMTWaveformSimConfig | 6 + .../PMTfittingparametersWhitespace.csv | 122 +++++++ .../PMTWaveformSim/PhaseIIADCHitFinderConfig | 11 + configfiles/PMTWaveformSim/README.md | 25 ++ configfiles/PMTWaveformSim/ToolChainConfig | 26 ++ configfiles/PMTWaveformSim/ToolsConfig | 4 + 18 files changed, 816 insertions(+), 370 deletions(-) create mode 100644 UserTools/PMTWaveformSim/PMTWaveformSim.cpp create mode 100644 UserTools/PMTWaveformSim/PMTWaveformSim.h create mode 100644 UserTools/PMTWaveformSim/README.md mode change 100755 => 100644 UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp mode change 100755 => 100644 UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.h create mode 100644 configfiles/PMTWaveformSim/DummyToolConfig create mode 100644 configfiles/PMTWaveformSim/LoadGeometryConfig create mode 100644 configfiles/PMTWaveformSim/LoadWCSimConfig create mode 100644 configfiles/PMTWaveformSim/PMTWaveformSimConfig create mode 100644 configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv create mode 100644 configfiles/PMTWaveformSim/PhaseIIADCHitFinderConfig create mode 100644 configfiles/PMTWaveformSim/README.md create mode 100644 configfiles/PMTWaveformSim/ToolChainConfig create mode 100644 configfiles/PMTWaveformSim/ToolsConfig diff --git a/DataModel/Waveform.h b/DataModel/Waveform.h index 28d1b71a3..f564a0d8b 100644 --- a/DataModel/Waveform.h +++ b/DataModel/Waveform.h @@ -55,4 +55,58 @@ class Waveform : public SerialisableObject{ } }; +template +class MCWaveform : public Waveform { + + friend class boost::serialization::access; + + public: + MCWaveform() : Waveform(), fParents(std::vector>{}) {} + MCWaveform(double tsin, std::vector samplesin, std::vector> theparents) : Waveform(tsin, samplesin), fParents(theparents) {} + virtual ~MCWaveform(){}; + + inline const std::vector>* GetParents() { return &fParents; } + inline const std::vector>& Parents() { return fParents; } + inline std::vector ParentsAtSample(int sample) { return fParents[sample]; } + + inline void SetParents(std::vector> parentsin){ fParents = parentsin; } + + inline Waveform GetBaseWaveform(){ return *this; } + + // Override the base print to add in parent info + bool Print() { + int verbose=0; + cout << "StartTime : " << this->fStartTime << endl; + cout << "NSamples : " << this->fSamples.size() << endl; + if(verbose){ + cout << "Samples (parents) : {"; + for(int samplei=0; samplei < this->fSamples.size(); samplei++){ + cout << this->fSamples.at(samplei); + + if (this->fParents.size() == this->fSamples.size()) { + cout << "("; + for (auto parent : this->ParentsAtSample(samplei)) { + cout << parent; + if ((samplei+1) != this->ParentsAtSample(samplei).size()) + cout << ", "; + else + cout << ")"; + } + } + + if ((samplei+1) != this ->fSamples.size()) cout << ", "; + } + cout<<"}"<> fParents; + +}; + #endif diff --git a/UserTools/Factory/Factory.cpp b/UserTools/Factory/Factory.cpp index 9a06a5520..805aad1be 100644 --- a/UserTools/Factory/Factory.cpp +++ b/UserTools/Factory/Factory.cpp @@ -1,172 +1,12 @@ #include "Factory.h" Tool* Factory(std::string tool) { -Tool* ret=0; + Tool* ret = 0; + if (tool=="PlotWaveforms") ret=new PlotWaveforms; + if (tool=="LoadGeometry") ret=new LoadGeometry; + if (tool=="LoadWCSim") ret=new LoadWCSim; + if (tool=="PMTWaveformSim") ret=new PMTWaveformSim; + if (tool=="PhaseIIADCHitFinder") ret=new PhaseIIADCHitFinder; -// if (tool=="Type") tool=new Type; -if (tool=="DummyTool") ret=new DummyTool; -if (tool=="ExampleGenerateData") ret=new ExampleGenerateData; -if (tool=="ExampleSaveStore") ret=new ExampleSaveStore; -if (tool=="ExampleSaveRoot") ret=new ExampleSaveRoot; -if (tool=="ExampleloadStore") ret=new ExampleloadStore; -if (tool=="ExamplePrintData") ret=new ExamplePrintData; -if (tool=="ExampleLoadRoot") ret=new ExampleLoadRoot; -if (tool=="LAPPDParseScope") ret=new LAPPDParseScope; -if (tool=="LAPPDParseACC") ret=new LAPPDParseACC; -if (tool=="LAPPDFindPeak") ret=new LAPPDFindPeak; -if (tool=="SaveANNIEEvent") ret=new SaveANNIEEvent; -if (tool=="LAPPDSim") ret=new LAPPDSim; -if (tool=="LoadWCSim") ret=new LoadWCSim; -if (tool=="FindMrdTracks") ret=new FindMrdTracks; -if (tool=="PrintANNIEEvent") ret=new PrintANNIEEvent; -if (tool=="GenerateHits") ret=new GenerateHits; -if (tool=="LAPPDcfd") ret=new LAPPDcfd; -if (tool=="LAPPDBaselineSubtract") ret=new LAPPDBaselineSubtract; -if (tool=="NeutronStudyReadSandbox") ret=new NeutronStudyReadSandbox; -if (tool=="NeutronStudyPMCS") ret=new NeutronStudyPMCS; -if (tool=="NeutronStudyWriteTree") ret=new NeutronStudyWriteTree; -if (tool=="BeamTimeAna") ret=new BeamTimeAna; -if (tool=="BeamTimeTreeMaker") ret=new BeamTimeTreeMaker; -if (tool=="BeamTimeTreeReader") ret=new BeamTimeTreeReader; -if (tool=="RawLoader") ret=new RawLoader; -if (tool=="LAPPDBaselineSubtract") ret=new LAPPDBaselineSubtract; -if (tool=="LAPPDSaveROOT") ret=new LAPPDSaveROOT; -if (tool=="LAPPDFilter") ret=new LAPPDFilter; -if (tool=="LAPPDIntegratePulse") ret=new LAPPDIntegratePulse; -if (tool=="ADCCalibrator") ret=new ADCCalibrator; -if (tool=="ADCHitFinder") ret=new ADCHitFinder; -if (tool=="BeamChecker") ret=new BeamChecker; -if (tool=="BeamFetcher") ret=new BeamFetcher; -if (tool=="FindTrackLengthInWater") ret=new FindTrackLengthInWater; -if (tool=="LoadANNIEEvent") ret=new LoadANNIEEvent; -if (tool=="PhaseITreeMaker") ret=new PhaseITreeMaker; -if (tool=="MrdPaddlePlot") ret=new MrdPaddlePlot; -if (tool=="LoadWCSimLAPPD") ret=new LoadWCSimLAPPD; -if (tool=="WCSimDemo") ret=new WCSimDemo; -if (tool=="DigitBuilder") ret=new DigitBuilder; -if (tool=="VtxSeedGenerator") ret=new VtxSeedGenerator; -if (tool=="VtxPointPositionFinder") ret=new VtxPointPositionFinder; -if (tool=="LAPPDlasertestHitFinder") ret=new LAPPDlasertestHitFinder; -if (tool=="RawLoadToRoot") ret=new RawLoadToRoot; -if (tool=="MRDPulseFinder") ret=new MRDPulseFinder; -if (tool=="LAPPDAnalysis") ret=new LAPPDAnalysis; -if (tool=="ExampleOverTool") ret=new ExampleOverTool; -if (tool=="PhaseIITreeMaker") ret=new PhaseIITreeMaker; -if (tool=="VertexGeometryCheck") ret=new VertexGeometryCheck; -if (tool=="LikelihoodFitterCheck") ret=new LikelihoodFitterCheck; -if (tool=="EventSelector") ret=new EventSelector; -if (tool=="SaveRecoEvent") ret=new SaveRecoEvent; -if (tool=="VtxExtendedVertexFinder") ret=new VtxExtendedVertexFinder; -if (tool=="VtxPointDirectionFinder") ret=new VtxPointDirectionFinder; -if (tool=="VtxPointVertexFinder") ret=new VtxPointVertexFinder; -if (tool=="LoadCCData") ret=new LoadCCData; -if (tool=="HitCleaner") ret=new HitCleaner; -if (tool=="HitResiduals") ret=new HitResiduals; -if (tool=="MonitorReceive") ret=new MonitorReceive; -if (tool=="MonitorSimReceive") ret=new MonitorSimReceive; -if (tool=="DigitBuilderDoE") ret=new DigitBuilderDoE; -if (tool=="EventSelectorDoE") ret=new EventSelectorDoE; -if (tool=="MonitorMRDTime") ret=new MonitorMRDTime; -if (tool=="MonitorMRDLive") ret=new MonitorMRDLive; -if (tool=="PulseSimulation") ret=new PulseSimulation; -if (tool=="PlotLAPPDTimesFromStore") ret=new PlotLAPPDTimesFromStore; -if (tool=="CheckDetectorCounts") ret=new CheckDetectorCounts; -if (tool=="MrdDistributions") ret=new MrdDistributions; -if (tool=="MCParticleProperties") ret=new MCParticleProperties; -if (tool=="DigitBuilderROOT") ret=new DigitBuilderROOT; -if (tool=="MrdEfficiency") ret=new MrdEfficiency; -if (tool=="EventDisplay") ret=new EventDisplay; -if (tool=="TankCalibrationDiffuser") ret=new TankCalibrationDiffuser; -if (tool=="TotalLightMap") ret=new TotalLightMap; -if (tool=="MrdDiscriminatorScan") ret=new MrdDiscriminatorScan; -if (tool=="MCRecoEventLoader") ret=new MCRecoEventLoader; -if (tool=="MonitorMRDEventDisplay") ret=new MonitorMRDEventDisplay; -if (tool=="LoadGeometry") ret=new LoadGeometry; -if (tool=="LoadRATPAC") ret=new LoadRATPAC; -if (tool=="TimeClustering") ret=new TimeClustering; -if (tool=="GracefulStop") ret=new GracefulStop; -if (tool=="PhaseIIADCHitFinder") ret=new PhaseIIADCHitFinder; -if (tool=="TrackCombiner") ret=new TrackCombiner; -if (tool=="SimulatedWaveformDemo") ret=new SimulatedWaveformDemo; -if (tool=="CNNImage") ret=new CNNImage; -if (tool=="MonitorTankTime") ret=new MonitorTankTime; -if (tool=="PhaseIIADCCalibrator") ret=new PhaseIIADCCalibrator; -if (tool=="MCHitToHitComparer") ret=new MCHitToHitComparer; -if (tool=="MCPropertiesToTree") ret=new MCPropertiesToTree; -if (tool=="CalcClassificationVars") ret=new CalcClassificationVars; -if (tool=="StoreClassificationVars") ret=new StoreClassificationVars; -if (tool=="LoadGenieEvent") ret=new LoadGenieEvent; -if (tool=="PrintGenieEvent") ret=new PrintGenieEvent; -if (tool=="PlotWaveforms") ret=new PlotWaveforms; -if (tool=="PMTDataDecoder") ret=new PMTDataDecoder; -if (tool=="ANNIEEventBuilder") ret=new ANNIEEventBuilder; -if (tool=="MRDDataDecoder") ret=new MRDDataDecoder; -if (tool=="PrintADCData") ret=new PrintADCData; -if (tool=="ClusterFinder") ret=new ClusterFinder; -if (tool=="PrintRecoEvent") ret=new PrintRecoEvent; -if (tool=="RunValidation") ret=new RunValidation; -if (tool=="AmBeRunStatistics") ret=new AmBeRunStatistics; -if (tool=="SimpleTankEnergyCalibrator") ret=new SimpleTankEnergyCalibrator; -if (tool=="BeamClusterPlots") ret=new BeamClusterPlots; -if (tool=="MrdPaddleEfficiencyPreparer") ret=new MrdPaddleEfficiencyPreparer; -if (tool=="MrdPaddleEfficiencyCalc") ret=new MrdPaddleEfficiencyCalc; -if (tool=="FMVEfficiency") ret=new FMVEfficiency; -if (tool=="LoadRawData") ret=new LoadRawData; -if (tool=="TriggerDataDecoder") ret=new TriggerDataDecoder; -if (tool=="ClusterClassifiers") ret=new ClusterClassifiers; -if (tool=="MRDLoopbackAnalysis") ret=new MRDLoopbackAnalysis; -if (tool=="VetoEfficiency") ret=new VetoEfficiency; -if (tool=="EnergyExtractor") ret=new EnergyExtractor; -if (tool=="MonitorTrigger") ret=new MonitorTrigger; -if (tool=="LAPPDPlotWaveForms") ret=new LAPPDPlotWaveForms; -if (tool=="LAPPDReorderData") ret=new LAPPDReorderData; -if (tool=="LAPPDStoreReorderData") ret=new LAPPDStoreReorderData; -if (tool=="LAPPDFindT0") ret=new LAPPDFindT0; -if (tool=="LAPPDStoreFindT0") ret=new LAPPDStoreFindT0; -if (tool=="LAPPDStoreReadIn") ret=new LAPPDStoreReadIn; -if (tool=="ClusterDummy") ret=new ClusterDummy; -if (tool=="LAPPDCluster") ret=new LAPPDCluster; -if (tool=="LAPPDMakePeds") ret=new LAPPDMakePeds; -if (tool=="EventClassification") ret=new EventClassification; -if (tool=="DataSummary") ret=new DataSummary; -if (tool=="MaxPEPlots") ret=new MaxPEPlots; -if (tool=="StoreDecodedTimestamps") ret=new StoreDecodedTimestamps; -if (tool=="PlotDecodedTimestamps") ret=new PlotDecodedTimestamps; -if (tool=="MonitorDAQ") ret=new MonitorDAQ; -if (tool=="MonitorSimReceiveLAPPD") ret=new MonitorSimReceiveLAPPD; -if (tool=="MonitorLAPPDSC") ret=new MonitorLAPPDSC; -if (tool=="MonitorLAPPDData") ret=new MonitorLAPPDData; -if (tool=="ParseDataMonitoring") ret=new ParseDataMonitoring; -if (tool=="LAPPDClusterTree") ret=new LAPPDClusterTree; -if (tool=="LAPPDPlotWaveForms2D") ret=new LAPPDPlotWaveForms2D; -if (tool=="LAPPDGausBaselineSubtraction") ret=new LAPPDGausBaselineSubtraction; -if (tool=="LAPPDASCIIReadIn") ret=new LAPPDASCIIReadIn; -if (tool=="BeamDecoder") ret=new BeamDecoder; -if (tool=="LoadRunInfo") ret=new LoadRunInfo; -if (tool=="ApplyMRDEff") ret=new ApplyMRDEff; -if (tool=="SimpleReconstruction") ret=new SimpleReconstruction; -if (tool=="LAPPDnnlsPeak") ret=new LAPPDnnlsPeak; -if (tool=="LAPPDLocateHit") ret=new LAPPDLocateHit; -if (tool=="LAPPDOtherSimp") ret=new LAPPDOtherSimp; -if (tool=="LAPPDTraceMax") ret=new LAPPDTraceMax; -if (tool=="saveLAPPDInfo") ret=new saveLAPPDInfo; -if (tool=="parseLAPPDData") ret=new parseLAPPDData; -if (tool=="checkLAPPDStatus") ret=new checkLAPPDStatus; -if (tool=="GetLAPPDEvents") ret=new GetLAPPDEvents; -if (tool=="LAPPDDataDecoder") ret=new LAPPDDataDecoder; -if (tool=="PythonScript") ret=new PythonScript; -if (tool=="MeanTimeCheck") ret=new MeanTimeCheck; -if (tool=="LoadReweightGenieEvent") ret=new LoadReweightGenieEvent; -if (tool=="ReweightFlux") ret=new ReweightFlux; -if (tool=="FilterLAPPDEvents") ret=new FilterLAPPDEvents; -if (tool=="VtxSeedFineGrid") ret=new VtxSeedFineGrid; -if (tool=="FilterEvents") ret=new FilterEvents; -if (tool=="Stage1DataBuilder") ret=new Stage1DataBuilder; -if (tool=="BeamFetcherV2") ret=new BeamFetcherV2; -if (tool=="FindNeutrons") ret=new FindNeutrons; -if (tool=="NeutronMultiplicity") ret=new NeutronMultiplicity; -if (tool=="PlotsTrackLengthAndEnergy") ret=new PlotsTrackLengthAndEnergy; -if (tool=="SaveConfigInfo") ret=new SaveConfigInfo; -if (tool=="ReadConfigInfo") ret=new ReadConfigInfo; -return ret; + return ret; } diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp new file mode 100644 index 000000000..9241bdfca --- /dev/null +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp @@ -0,0 +1,343 @@ +#include +#include +#include + +// ANNIE includes +#include "ANNIEconstants.h" +#include "PMTWaveformSim.h" + +// ROOT includes +#include "TGraph.h" + +PMTWaveformSim::PMTWaveformSim():Tool(){} + + +bool PMTWaveformSim::Initialise(std::string configfile, DataModel &data) +{ + + /////////////////// Useful header /////////////////////// + if(configfile!="") m_variables.Initialise(configfile); // loading config file + //m_variables.Print(); + + m_data= &data; //assigning transient data pointer + ///////////////////////////////////////////////////////////////// + + // get config variables + bool gotVerbosity = m_variables.Get("verbosity", verbosity); + if (!gotVerbosity) verbosity = 1; + + bool gotPMTParamFile = m_variables.Get("PMTParameterFile", fPMTParameterFile); + if (!gotPMTParamFile) { + logmessage = "PMTWaveformSim: No PMTParameterFile specified! Aborting!"; + Log(logmessage, v_error, verbosity); + return false; + } + + if (!LoadPMTParameters()) + return false; + + bool gotPrewindow = m_variables.Get("Prewindow", fPrewindow); + if (!gotPrewindow) { + logmessage = "PMTWaveformSim: Prewindow not defined. Using default of 10."; + Log(logmessage, v_warning, verbosity); + fPrewindow = 10; + } + + bool gotReadoutWindow = m_variables.Get("ReadoutWindow", fReadoutWindow); + if (!gotReadoutWindow) { + logmessage = "PMTWaveformSim: ReadoutWindow not defined. Using default of 35."; + Log(logmessage, v_warning, verbosity); + fReadoutWindow = 35; + } + + bool gotT0Offset = m_variables.Get("T0Offset", fT0Offset); + if (!gotT0Offset) { + logmessage = "PMTWaveformSim: T0Offset not defined. Using default of 7."; + Log(logmessage, v_warning, verbosity); + fT0Offset = 7; + } + + bool gotDebug = m_variables.Get("MakeDebugFile", fDebug); + if (!gotDebug) fDebug = 0; + if (fDebug) + fOutFile = new TFile("PMTWaveforms.root", "RECREATE"); + + bool gotGeometry = m_data->Stores.at("ANNIEEvent")->Header->Get("AnnieGeometry", fGeo); + if(!gotGeometry){ + logmessage = "PMTWaveformSim: Error retrieving Geometry from ANNIEEvent! Aborting!"; + Log(logmessage, v_error, verbosity); + return false; + } + + // Could set a seed here for repeatability + // Though would probably want to seed it in the Execute function based on the run/part/event numbers + fRandom = new TRandom3(); + + return true; +} + +//------------------------------------------------------------------------------ +bool PMTWaveformSim::Execute() +{ + + if (!LoadFromStores()) + return false; + + // The container for the data that we'll put into the ANNIEEvent + std::map> > RawADCDataMC; + std::map> > CalADCDataMC; + + for (auto mcHitsIt : *fMCHits) { // Loop over the hit PMTs + int PMTID = mcHitsIt.first; + + std::cout << "PMT: " << PMTID << std::endl; + + std::vector mcHits = mcHitsIt.second; + + // Generate waveform samples from the MC hits + // samples from hits that are close in time will be added together + // key is hit time in clock ticks, value is amplitude + std::map sample_map; + std::map> parent_map; // use a set so each parent is recorded only once + + for (const auto& mcHit : mcHits) { // Loop through each MCHit in the vector + // convert PMT hit time to clock ticks and "digitize" by converting to an int + // skip negative hit times, what does that even mean if we're not using the smeared digit time? + if (mcHit.GetTime() < 0) continue; + + // MCHit time is in ns, but we're going to sample in clock ticks + uint16_t hit_t0 = uint16_t(mcHit.GetTime() / NS_PER_ADC_SAMPLE); + double hit_charge = mcHit.GetCharge(); + + std::cout << mcHit.GetTime() << ", " << hit_charge << std::endl; + + // Set the readout window but don't allow negative times + uint16_t start_clocktick = (hit_t0 > fPrewindow)? hit_t0 - fPrewindow : 0; + uint16_t end_clocktick = start_clocktick + fReadoutWindow; + + // loop over clock ticks + for (uint16_t clocktick = start_clocktick; clocktick <= end_clocktick; clocktick += 1) { + uint16_t sample = CustomLogNormalPulse(hit_t0, clocktick, hit_charge, PMTID); + + // check if this hit time has been recorded + // either set it or add to it + if (sample_map.find(clocktick) == sample_map.end()) { + sample_map[clocktick] = sample; + parent_map[clocktick] = std::set(mcHit.GetParents()->begin(), mcHit.GetParents()->end()); + } else { + sample_map[clocktick] += sample; + + // Put the parents into the set + for (uint idx = 0; idx < mcHit.GetParents()->size(); ++idx) + parent_map[clocktick].insert(mcHit.GetParents()->at(idx)); + } + }// end loop over clock ticks + }// end loop over mcHits + + // Set the noise envelope and baseline for this PMT + // The noise std dev appears to be normally distributed around 1 with sigma 0.25 + double noiseSigma = fRandom->Gaus(1, 0.25); + int basline = fRandom->Uniform(300, 350); + + // convert the sample map into a vector of Waveforms and put them into the container + std::vector> rawWaveforms; + std::vector> calWaveforms; + ConvertMapToWaveforms(sample_map, parent_map, rawWaveforms, calWaveforms, noiseSigma, basline); + + RawADCDataMC.emplace(PMTID, rawWaveforms); + CalADCDataMC.emplace(PMTID, calWaveforms); + }// end loop over PMTs + + // Publish the waveforms to the ANNIEEvent store + m_data->Stores.at("ANNIEEvent")->Set("RawADCDataMC", RawADCDataMC); + m_data->Stores.at("ANNIEEvent")->Set("CalibratedADCData", CalADCDataMC); + + if (fDebug) + FillDebugGraphs(RawADCDataMC); + + ++fEvtNum; + return true; +} + +//------------------------------------------------------------------------------ +bool PMTWaveformSim::Finalise() +{ + if (fDebug) + fOutFile->Close(); + + return true; +} + +//------------------------------------------------------------------------------ +bool PMTWaveformSim::LoadPMTParameters() +{ + + std::ifstream infile(fPMTParameterFile); + if (!infile.is_open()) { + logmessage = "PMTWaveformSim: Error opening CSV file: "; + logmessage += fPMTParameterFile + "!"; + Log(logmessage, v_error, verbosity); + return false; + } + + int pmtid; + double p0, p1, p2; + std::string line; + std::getline(infile, line); // Skipping the header line + + while (std::getline(infile, line)) { + if (infile.fail()) { + logmessage = "PMTWaveformSim: Error reading from CSV file: "; + logmessage += fPMTParameterFile + "!"; + Log(logmessage, v_error, verbosity); + return false; + } + + // Skip any commented lines + if(line.find("#")!=std::string::npos) continue; + + // Turn the line into a stringstream to extract the values + std::stringstream ss(line); + ss >> pmtid >> p0 >> p1 >> p2; + + pmtParameters[pmtid] = std::make_tuple(p0, p1, p2); + + logmessage = "PMTWaveformSim: Loaded parameters for PMTID " + std::to_string(pmtid) + ": "; + logmessage += "p0 = " + std::to_string(p0); + logmessage += "p1 = " + std::to_string(p1); + logmessage += "p2 = " + std::to_string(p2); + Log(logmessage, v_message, verbosity); + } + + infile.close(); + + return true; +} + +//------------------------------------------------------------------------------ +uint16_t PMTWaveformSim::CustomLogNormalPulse(uint16_t hit_t0, uint16_t clocktick, double hit_charge, int PMTID) +{ + double p0, p1, p2; + if (pmtParameters.find(PMTID) != pmtParameters.end()) { + std::tie(p0, p1, p2) = pmtParameters[PMTID]; + } else { + logmessage = "PMTWaveformSim: PMTParameters not found for " + std::to_string(PMTID); + logmessage += ", using defaults: p0 = 17.49, p1 = 3.107, p2 = 0.1492"; + Log(logmessage, v_warning, verbosity); + + p0 = 17.49; + p1 = 3.107; + p2 = 0.1492; + } + + // The fit was performed in time units of ns, but we want to grab samples in clock ticks + double x = (double(clocktick) + fT0Offset - hit_t0) * NS_PER_ADC_SAMPLE; + // std::cout << " " << clocktick << ", " << x << std::endl; + + double numerator = pow(log(x) - p1, 2); + double denom = (2 * pow(p2, 2)); + double amplitude = p0 * exp(-numerator / denom) * hit_charge; + + // Clip at 4095 and digitize to an integer + return uint16_t((amplitude > 4095) ? 4095 : amplitude); +} + +//------------------------------------------------------------------------------ +void PMTWaveformSim::ConvertMapToWaveforms(const std::map &sample_map, + const std::map> & parent_map, + std::vector> &rawWaveforms, + std::vector> &calWaveforms, + double noiseSigma, int baseline) +{ + // All MC has extended readout + std::vector rawSamples; + std::vector calSamples; + std::vector> parents; + for (uint16_t tick = 0; tick < 34993; ++tick) { + // Generate noise for each sample based on the std dev of the noise envelope + double noise = fRandom->Gaus(0, noiseSigma); + + int sample = std::round(noise + baseline); + + // look for this tick in the sample map and add it + if (sample_map.find(tick) != sample_map.end()) + sample += sample_map.at(tick); + + + rawSamples.push_back((sample > 4095) ? 4095 : sample); + calSamples.push_back((rawSamples.back() - baseline) * ADC_TO_VOLT); + + std::vector innerParents; + if (parent_map.find(tick) != parent_map.end()) + std::copy(parent_map.at(tick).begin(), parent_map.at(tick).begin(), + std::back_inserter(innerParents)); + else + innerParents.push_back(-5); + + parents.push_back(innerParents); + } + + // The start time in data is a timestamp. Don't have that for MC so just set to 0. + rawWaveforms.emplace_back(0, rawSamples, parents); + calWaveforms.emplace_back(0, calSamples, baseline, noiseSigma); +} + +//------------------------------------------------------------------------------ +bool PMTWaveformSim::LoadFromStores() +{ + bool goodAnnieEvent = m_data->Stores.count("ANNIEEvent"); + if (!goodAnnieEvent) { + logmessage = "PMTWaveformSim: no ANNIEEvent store!"; + Log(logmessage, v_error, verbosity); + return false; + } + + bool goodMCHits = m_data->Stores.at("ANNIEEvent")->Get("MCHits", fMCHits); + if (!goodMCHits) { + logmessage = "PMTWaveformSim: no MCHits in the ANNIEEvent!"; + Log(logmessage, v_error, verbosity); + return false; + } + + return true; + +} + +//------------------------------------------------------------------------------ +void PMTWaveformSim::FillDebugGraphs(const std::map> > &RawADCDataMC) +{ + for (auto itpair : RawADCDataMC) { + std::string chanString = std::to_string(itpair.first); + + // Get/make the directory for this PMT + TDirectory* dir = fOutFile->GetDirectory(chanString.c_str()); + if (!dir) + dir = fOutFile->mkdir(chanString.c_str()); + + // Hop into that directory and save the graph + dir->cd(); + + // loop over waveforms and make graphs + for (uint wfIdx = 0; wfIdx < itpair.second.size(); ++wfIdx) { + auto waveform = itpair.second.at(wfIdx); + + std::string grName = ("wf_" + std::to_string(fEvtNum) + "_" + std::to_string(wfIdx)); + + // Make the graph + std::vector samples = waveform.Samples(); + TGraph* grTemp = new TGraph(); + double sampleX = waveform.GetStartTime() / NS_PER_ADC_SAMPLE; + for(auto sample : samples) { + grTemp->AddPoint(sampleX, sample); + ++sampleX; + } + + grTemp->Write(grName.c_str()); + }// end loop over waveforms + }// end loop over PMTs +} + + + + + diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.h b/UserTools/PMTWaveformSim/PMTWaveformSim.h new file mode 100644 index 000000000..c46d8bf75 --- /dev/null +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.h @@ -0,0 +1,75 @@ +#ifndef PMTWaveformSim_H +#define PMTWaveformSim_H + +#include +#include + +// ANNIE includes +#include "Tool.h" +#include "CalibratedADCWaveform.h" +#include "Waveform.h" + +// ROOT includes +#include "TFile.h" +#include "TRandom3.h" + + +/** + * \class PMTWaveformSim + * + * This is a blank template for a Tool used by the script to generate a new custom tool. Please fill out the description and author information. +* +* $Author: D. Ajana $ +* $Date: 2024/11/05 10:44:00 $ +* Contact: dja23@fsu.edu +*/ +class PMTWaveformSim: public Tool { + + + public: + + PMTWaveformSim(); ///< Simple constructor + bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools. + bool Execute(); ///< Execute function used to perform Tool purpose. + bool Finalise(); ///< Finalise function used to clean up resources. + bool LoadFromStores(); + + bool LoadPMTParameters(); + uint16_t CustomLogNormalPulse(uint16_t hit_t0, uint16_t t0_clocktick, double hit_charge, int PMTID); + void ConvertMapToWaveforms(const std::map &sample_map, + const std::map> & parent_map, + std::vector> &rawWaveforms, + std::vector> &calWaveforms, + double noiseSigma, int baseline); + + void FillDebugGraphs(const std::map> > &RawADCDataMC); + + private: + + // To load from the ANNIEEvent + std::map> *fMCHits = nullptr; + std::map> pmtParameters; + Geometry *fGeo = nullptr; + + // Config variables + uint16_t fPrewindow; + uint16_t fReadoutWindow; + uint16_t fT0Offset; + std::string fPMTParameterFile; + + TRandom3 *fRandom; + + bool fDebug; + TFile *fOutFile; + int fEvtNum = 0; + + int verbosity; + int v_error=0; + int v_warning=1; + int v_message=2; + int v_debug=3; + std::string logmessage; +}; + + +#endif diff --git a/UserTools/PMTWaveformSim/README.md b/UserTools/PMTWaveformSim/README.md new file mode 100644 index 000000000..4142083da --- /dev/null +++ b/UserTools/PMTWaveformSim/README.md @@ -0,0 +1,20 @@ +# PMTWaveformSim + +PMTWaveformSim + +## Data + +Describe any data formats PMTWaveformSim creates, destroys, changes, or analyzes. E.G. + +**RawLAPPDData** `map>>` +* Takes this data from the `ANNIEEvent` store and finds the number of peaks + + +## Configuration + +Describe any configuration variables for PMTWaveformSim. + +``` +param1 value1 +param2 value2 +``` diff --git a/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp b/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp index a85341ac1..82c5d45fc 100644 --- a/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp +++ b/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp @@ -723,20 +723,27 @@ PhaseIIADCCalibrator::make_calibrated_waveforms_ze3ra_multi( // come from the same readout for the same channel) std::vector< CalibratedADCWaveform > calibrated_waveforms; for (const auto& raw_waveform : raw_waveforms) { - std::vector baselines; + std::vector baselines; + std::vector sigma_baselines; std::vector RepresentationRegion; - double first_baseline=0; - double first_sigma_baseline; + bool isFirst = true; + double first_baseline, first_sigma_baseline; double baseline, sigma_baseline; const size_t nsamples = raw_waveform.Samples().size(); for(size_t starting_sample = 0; starting_sample < nsamples; starting_sample += baseline_rep_samples){ double baseline, sigma_baseline; ze3ra_baseline(raw_waveform, baseline, sigma_baseline, num_baseline_samples,starting_sample); - if(sigma_baseline4) std::cout << "BASELINE UNCERTAINTY BEYOND SET THRESHOLD. IGNORING SAMPLE" << std::endl; } } @@ -746,6 +753,7 @@ PhaseIIADCCalibrator::make_calibrated_waveforms_ze3ra_multi( if(verbosity>4) std::cout << "NO BASLINE FOUND WITHIN TOLERANCE. USING FIRST AS BEST ESTIMATE" << std::endl; RepresentationRegion.push_back(baseline_rep_samples); baselines.push_back(first_baseline); + sigma_baselines.push_back(first_sigma_baseline); } std::vector cal_data; const std::vector& raw_data = raw_waveform.Samples(); @@ -762,10 +770,22 @@ PhaseIIADCCalibrator::make_calibrated_waveforms_ze3ra_multi( } } } - double bl_estimates_mean, bl_estimates_var; - ComputeMeanAndVariance(baselines, bl_estimates_mean, bl_estimates_var); + double bl_estimates_mean = 0; + double bl_estimates_sigma = 0; + + // When averaging multiple means you need to average the variances as well + // We don't want to use the variance of the mean of all the baselines. + // If there is only one baseline then the variance is 0 if you do that... + for (size_t idx = 0; idx < baselines.size(); ++idx) { + bl_estimates_mean += baselines[idx]; + bl_estimates_sigma += pow(sigma_baselines[idx],2); + } + bl_estimates_sigma = sqrt(bl_estimates_sigma / double(baselines.size())); + + // ComputeMeanAndVariance(baselines, bl_estimates_mean, bl_estimates_var, std::numeric_limits::max(), 0, 7); + calibrated_waveforms.emplace_back(raw_waveform.GetStartTime(), - cal_data, bl_estimates_mean, bl_estimates_var); + cal_data, bl_estimates_mean, bl_estimates_sigma); } return calibrated_waveforms; } diff --git a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp old mode 100755 new mode 100644 index 8be738345..cb05bd7d2 --- a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp +++ b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp @@ -36,6 +36,7 @@ bool PhaseIIADCHitFinder::Initialise(std::string config_filename, DataModel& dat m_variables.Get("PulseWindowEnd", pulse_window_end_shift); m_variables.Get("WindowIntegrationDB", adc_window_db); m_variables.Get("EventBuilding",eventbuilding_mode); + m_variables.Get("MCWaveforms",mc_waveforms); if ((pulse_window_start_shift > 0) || (pulse_window_end_shift) < 0){ Log("PhaseIIADCHitFinder Tool: WARNING... trigger threshold crossing will not be inside pulse window. Threshold" @@ -65,6 +66,11 @@ bool PhaseIIADCHitFinder::Initialise(std::string config_filename, DataModel& dat eventbuilding_mode = false; } + if (mc_waveforms && (eventbuilding_mode || use_led_waveforms)) { + Log("PhaseIIADCCalibrator: Cannot use MCWaveforms in EventBuilding mode or while using LED waveforms. Aborting!", v_error, verbosity); + return false; + } + //Set in CStore for tools to know and log this later m_data->CStore.Set("ADCThreshold",default_adc_threshold); @@ -129,22 +135,50 @@ bool PhaseIIADCHitFinder::Execute() { } else { got_raw_data = annie_event->Get("RawADCData", raw_waveform_map); got_rawaux_data = annie_event->Get("RawADCAuxData", raw_aux_waveform_map); + + if (mc_waveforms) { + raw_waveform_map.clear(); + + std::map > > mc_waveform_map; + bool got_raw_mc = annie_event->Get("RawADCDataMC", mc_waveform_map); + if (!got_raw_mc) { + Log("Error: The PhaseIIADCHitFinder tool could not find the RawADCDataMC entry", v_error, + verbosity); + return false; + } + + // Slice off the MC info to get a map of Waveforms + // Need to grab the parents for use in BackTracker later + for (auto itpair : mc_waveform_map) { + for (auto mc_waveform : itpair.second) + raw_waveform_map[itpair.first].push_back(mc_waveform.GetBaseWaveform()); + } + + if (raw_waveform_map.size() == mc_waveform_map.size()) + got_raw_data = true; + else { + Log("Error: The PhaseIIADCHitFinder tool could not extract Waveforms from the MCWaveforms.", v_error, + verbosity); + return false; + } + }// end if mc_waveforms } + // Check for problems if ( !got_raw_data ) { Log("Error: The PhaseIIADCHitFinder tool could not find the RawADCData entry", v_error, verbosity); return false; } - if ( !got_rawaux_data ) { + if ( !got_rawaux_data && !(use_led_waveforms || mc_waveforms)) { Log("Error: The PhaseIIADCHitFinder tool could not find the RawADCAuxData entry", v_error, - verbosity); + verbosity); return false; } else if ( raw_waveform_map.empty() ) { - Log("Error: The PhaseIIADCHitFinder tool found an empty RawADCData entry", v_error, - verbosity); - return false; + Log("Error: The PhaseIIADCHitFinder tool found an empty RawADCData entry. Skipping.", v_error, + verbosity); + return true; } // Load the maps containing the ADC calibrated waveform data @@ -165,15 +199,15 @@ bool PhaseIIADCHitFinder::Execute() { " entry", v_error, verbosity); return false; } - if ( !got_calibratedaux_data ) { + if ( !got_calibratedaux_data && !(use_led_waveforms || mc_waveforms)) { Log("Error: The PhaseIIADCHitFinder tool could not find the CalibratedADCAuxData" " entry", v_error, verbosity); return false; } else if ( calibrated_waveform_map.empty() ) { - Log("Error: The PhaseIIADCHitFinder tool found an empty CalibratedADCData entry", + Log("Error: The PhaseIIADCHitFinder tool found an empty CalibratedADCData entry. Skipping", v_error, verbosity); - return false; + return true; } //Find pulses in the raw detector data @@ -729,6 +763,8 @@ std::vector PhaseIIADCHitFinder::find_pulses_bywindow( } +// ****************************************************************** +// "PulseFindingApproach" default for EventBuilding std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( const Waveform& raw_minibuffer_data, const CalibratedADCWaveform& calibrated_minibuffer_data, @@ -824,8 +860,8 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( if(it != ChannelKeyToTimingOffsetMap.end()){ //Timing offset is available timing_offset = ChannelKeyToTimingOffsetMap.at(channel_key); } else { - if(verbosity>2){ - std::cout << "Didn't find Timing offset for channel " << channel_key << std::endl; + if(verbosity>v_error){ + std::cout << "PhaseIIADCHitFinder: Didn't find Timing offset for channel... setting this channel's offset to 0ns" << channel_key << std::endl; } } @@ -837,7 +873,10 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( calibrated_minibuffer_data.GetSigmaBaseline(), raw_area, max_ADC, calibrated_amplitude, charge); } - + + + // ****************************************************************** + // "PulseWindowType" default for EventBuilding // Peak windows are defined only by crossing and un-crossing of ADC threshold } else if(pulse_window_type == "dynamic"){ size_t pulse_start_sample = BOGUS_INT; @@ -898,8 +937,8 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( if(it != ChannelKeyToTimingOffsetMap.end()){ //Timing offset is available timing_offset = ChannelKeyToTimingOffsetMap.at(channel_key); } else { - if(verbosity>2){ - std::cout << "Didn't find Timing offset for channel " << channel_key << std::endl; + if(verbosity>v_error){ + std::cout << "PhaseIIADCHitFinder: Didn't find Timing offset for channel... setting this channel's offset to 0ns" << channel_key << std::endl; } } @@ -932,20 +971,23 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( } } - if(verbosity>4) std::cout << "Hit time [ns] " << hit_time * NS_PER_ADC_SAMPLE << std::endl; + if(verbosity>v_debug) { + + std::cout << "Hit time [ns] " << hit_time * NS_PER_ADC_SAMPLE << std::endl; - if (hit_time < 0.0) { - // If for some reason the interpolation finds a negative time value (if the pulse is extremely early in the buffer), - // default to the peak time - std::cout << "Hit time is negative! Defaulting to peak time" << std::endl; - hit_time = peak_sample; + if (hit_time < 0.0) { + // If for some reason the interpolation finds a negative time value (if the pulse is extremely early in the buffer), + // default to the peak time (maximum ADC value of the pulse) + std::cout << "Hit time is negative! Defaulting to peak time" << std::endl; + hit_time = peak_sample; + } } // Store the freshly made pulse in the vector of found pulses pulses.emplace_back(channel_key, ( pulse_start_sample * NS_PER_ADC_SAMPLE )-timing_offset, - (peak_sample * NS_PER_ADC_SAMPLE)-timing_offset, + (hit_time * NS_PER_ADC_SAMPLE)-timing_offset, // interpolated hit time calibrated_minibuffer_data.GetBaseline(), calibrated_minibuffer_data.GetSigmaBaseline(), raw_area, max_ADC, calibrated_amplitude, charge); diff --git a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.h b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.h old mode 100755 new mode 100644 index b3471d58e..eae2a98e8 --- a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.h +++ b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.h @@ -58,7 +58,8 @@ class PhaseIIADCHitFinder : public Tool { int pulse_window_end_shift; std::map channel_threshold_map; std::map>> channel_window_map; - bool eventbuilding_mode; + bool eventbuilding_mode; + bool mc_waveforms; std::map ChannelKeyToTimingOffsetMap; diff --git a/UserTools/Unity.h b/UserTools/Unity.h index 5f1e11fb7..ea6856348 100644 --- a/UserTools/Unity.h +++ b/UserTools/Unity.h @@ -1,178 +1,7 @@ -#include "DummyTool.h" -#include "ExampleGenerateData.h" -#include "ExampleSaveStore.h" -#include "ExampleSaveRoot.h" -#include "ExampleloadStore.h" -#include "ExamplePrintData.h" -#include "ExampleLoadRoot.h" -#include "LAPPDBaselineSubtract.h" -#include "LAPPDcfd.h" -#include "TSplineFit.h" -#include "TPoly3.h" -#include "TZigZag.h" -#include "TOnePadDisplay.h" -#include "TBandedLE.h" -#include "LAPPDFilter.h" -#include "LAPPDFindPeak.h" -#include "LAPPDIntegratePulse.h" -#include "LAPPDParseACC.h" -#include "LAPPDParseScope.h" -#include "LAPPDSaveROOT.h" -#include "SaveANNIEEvent.h" -#include "LAPPDSim.h" -#include "LAPPDresponse.h" -#include "BeamTimeAna.h" -#include "BeamTimeTreeMaker.h" -#include "BeamTimeTreeReader.h" -#include "FindMrdTracks.h" -#include "LoadWCSim.h" -#include "wcsimT.h" -#include "PrintANNIEEvent.h" -#include "GenerateHits.h" -#include "NeutronStudyReadSandbox.h" -#include "NeutronStudyPMCS.h" -#include "NeutronStudyWriteTree.h" -#include "RawLoader.h" -#include "HeftyTreeReader.h" #include "Unity_recoANNIE.h" -#include "ADCCalibrator.h" -#include "ADCHitFinder.h" -#include "BeamChecker.h" -#include "IFBeamDBInterface.h" -#include "BeamFetcher.h" -#include "FindTrackLengthInWater.h" -#include "LoadANNIEEvent.h" -#include "PhaseITreeMaker.h" -#include "MrdPaddlePlot.h" -#include "LoadWCSimLAPPD.h" -#include "LAPPDTree.h" -#include "WCSimDemo.h" -#include "DigitBuilder.h" -#include "VtxSeedGenerator.h" -#include "VtxPointPositionFinder.h" -#include "LAPPDlasertestHitFinder.h" -#include "RawLoadToRoot.h" -#include "MRDPulseFinder.h" -#include "LAPPDAnalysis.h" -#include "ExampleOverTool.h" -#include "PhaseIITreeMaker.h" -#include "VertexGeometryCheck.h" -#include "LikelihoodFitterCheck.h" -#include "EventSelector.h" -#include "SaveRecoEvent.h" -#include "VtxExtendedVertexFinder.h" -#include "VtxPointDirectionFinder.h" -#include "VtxPointVertexFinder.h" -#include "LoadCCData.h" -#include "MRDTree.h" -#include "PMTData.h" -//#include "LoadCCData/RunInformation.h" -//#include "LoadCCData/TrigData.h" -#include "HitCleaner.h" -#include "HitResiduals.h" -#include "MonitorReceive.h" -#include "MonitorSimReceive.h" -#include "EventSelectorDoE.h" -#include "DigitBuilderDoE.h" -#include "MonitorMRDTime.h" -#include "MonitorMRDLive.h" -#include "PulseSimulation.h" -#include "PlotLAPPDTimesFromStore.h" -#include "CheckDetectorCounts.h" -#include "MrdDistributions.h" -#include "MCParticleProperties.h" -#include "DigitBuilderROOT.h" -#include "MrdEfficiency.h" -#include "EventDisplay.h" -#include "TankCalibrationDiffuser.h" -#include "TotalLightMap.h" -#include "MrdDiscriminatorScan.h" -#include "MCRecoEventLoader.h" -#include "MonitorMRDEventDisplay.h" +#include "PlotWaveforms.h" +#include "Factory.h" #include "LoadGeometry.h" -#include "LoadRATPAC.h" -#include "TimeClustering.h" -#include "GracefulStop.h" +#include "LoadWCSim.h" +#include "PMTWaveformSim.h" #include "PhaseIIADCHitFinder.h" -#include "TrackCombiner.h" -#include "SimulatedWaveformDemo.h" -#include "CNNImage.h" -#include "MonitorTankTime.h" -#include "PhaseIIADCCalibrator.h" -#include "MCHitToHitComparer.h" -#include "MCPropertiesToTree.h" -#include "CalcClassificationVars.h" -#include "StoreClassificationVars.h" -#include "LoadGenieEvent.h" -#include "PrintGenieEvent.h" -#include "PlotWaveforms.h" -#include "PMTDataDecoder.h" -#include "ANNIEEventBuilder.h" -#include "MRDDataDecoder.h" -#include "PrintADCData.h" -#include "ClusterFinder.h" -#include "RunValidation.h" -#include "AmBeRunStatistics.h" -#include "SimpleTankEnergyCalibrator.h" -#include "BeamClusterPlots.h" -#include "PrintRecoEvent.h" -#include "MrdPaddleEfficiencyPreparer.h" -#include "MrdPaddleEfficiencyCalc.h" -#include "FMVEfficiency.h" -#include "LoadRawData.h" -#include "TriggerDataDecoder.h" -#include "ClusterClassifiers.h" -#include "MRDLoopbackAnalysis.h" -#include "VetoEfficiency.h" -#include "EnergyExtractor.h" -#include "MonitorTrigger.h" -#include "EventClassification.h" -#include "DataSummary.h" -#include "MaxPEPlots.h" -#include "StoreDecodedTimestamps.h" -#include "PlotDecodedTimestamps.h" -#include "MonitorDAQ.h" -#include "MonitorSimReceiveLAPPD.h" -#include "MonitorLAPPDSC.h" -#include "MonitorLAPPDData.h" -#include "ParseDataMonitoring.h" -#include "LAPPDReorderData.h" -#include "LAPPDStoreReorderData.h" -#include "LAPPDStoreReadIn.h" -#include "LAPPDPlotWaveForms.h" -#include "LAPPDFindT0.h" -#include "LAPPDStoreFindT0.h" -#include "ClusterDummy.h" -#include "LAPPDMakePeds.h" -#include "LAPPDCluster.h" -#include "LAPPDClusterTree.h" -#include "LAPPDPlotWaveForms2D.h" -#include "LAPPDGausBaselineSubtraction.h" -#include "LAPPDASCIIReadIn.h" -#include "BeamDecoder.h" -#include "LoadRunInfo.h" -#include "ApplyMRDEff.h" -#include "SimpleReconstruction.h" -#include "LAPPDnnlsPeak.h" -#include "LAPPDLocateHit.h" -#include "LAPPDOtherSimp.h" -#include "LAPPDTraceMax.h" -#include "saveLAPPDInfo.h" -#include "parseLAPPDData.h" -#include "checkLAPPDStatus.h" -#include "GetLAPPDEvents.h" -#include "LAPPDDataDecoder.h" -#include "PythonScript.h" -#include "MeanTimeCheck.h" -#include "LoadReweightGenieEvent.h" -#include "ReweightFlux.h" -#include "FilterLAPPDEvents.h" -#include "VtxSeedFineGrid.h" -#include "FilterEvents.h" -#include "Stage1DataBuilder.h" -#include "BeamFetcherV2.h" -#include "FindNeutrons.h" -#include "NeutronMultiplicity.h" -#include "PlotsTrackLengthAndEnergy.h" -#include "SaveConfigInfo.h" -#include "ReadConfigInfo.h" diff --git a/configfiles/PMTWaveformSim/DummyToolConfig b/configfiles/PMTWaveformSim/DummyToolConfig new file mode 100644 index 000000000..95cad88d2 --- /dev/null +++ b/configfiles/PMTWaveformSim/DummyToolConfig @@ -0,0 +1,3 @@ +# Dummy config file + +verbose 2 \ No newline at end of file diff --git a/configfiles/PMTWaveformSim/LoadGeometryConfig b/configfiles/PMTWaveformSim/LoadGeometryConfig new file mode 100644 index 000000000..b52737aa4 --- /dev/null +++ b/configfiles/PMTWaveformSim/LoadGeometryConfig @@ -0,0 +1,9 @@ +verbosity 0 +LAPPDChannelCount 60 +FACCMRDGeoFile ./configfiles/LoadGeometry/FullMRDGeometry.csv +DetectorGeoFile ./configfiles/LoadGeometry/DetectorGeometrySpecs.csv +LAPPDGeoFile ./configfiles/LoadGeometry/LAPPDGeometry.csv +TankPMTGeoFile ./configfiles/LoadGeometry/FullTankPMTGeometry.csv +TankPMTGainFile ./configfiles/LoadGeometry/ChannelSPEGains_BeamRun20192020.csv +AuxiliaryChannelFile ./configfiles/LoadGeometry/AuxChannels.csv +TankPMTTimingOffsetFile ./configfiles/LoadGeometry/TankPMTTimingOffsets.csv \ No newline at end of file diff --git a/configfiles/PMTWaveformSim/LoadWCSimConfig b/configfiles/PMTWaveformSim/LoadWCSimConfig new file mode 100644 index 000000000..5f6f6b0a9 --- /dev/null +++ b/configfiles/PMTWaveformSim/LoadWCSimConfig @@ -0,0 +1,16 @@ +#LoadWCSim Config File +# all variables retrieved with m_variables.Get() must be defined here! +verbose 0 + +InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standard/tank/pmt/wcsim_0.18.10.root + +WCSimVersion 3 ## should reflect the WCSim version of the files being loaded +HistoricTriggeroffset 0 ## time offset of digits relative to the trigger +UseDigitSmearedTime 0 ## whether to use smeared digit time (T), or true time of first photon (F) +LappdNumStrips 56 ## num channels to construct from each LAPPD +LappdStripLength 100 ## relative x position of each LAPPD strip, for dual-sided readout [mm] +LappdStripSeparation 10 ## stripline separation, for calculating relative y position of each LAPPD strip [mm] +PMTMask ./configfiles/LoadWCSim/DeadPMTIDs_p2v7.txt ## Which PMTs should be masked out? / are dead? +ChankeyToPMTIDMap ./configfiles/LoadWCSim/Chankey_WCSimID_v7.txt +ChankeyToMRDIDMap ./configfiles/LoadWCSim/MRD_Chankey_WCSimID.dat +ChankeyToFMVIDMap ./configfiles/LoadWCSim/FMV_Chankey_WCSimID.dat diff --git a/configfiles/PMTWaveformSim/PMTWaveformSimConfig b/configfiles/PMTWaveformSim/PMTWaveformSimConfig new file mode 100644 index 000000000..a271dfa82 --- /dev/null +++ b/configfiles/PMTWaveformSim/PMTWaveformSimConfig @@ -0,0 +1,6 @@ +verbosity 0 +PMTParameterFile configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv +Prewindow 10 +ReadoutWindow 35 +T0Offset 7 +MakeDebugFile 0 diff --git a/configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv b/configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv new file mode 100644 index 000000000..5d42d4eda --- /dev/null +++ b/configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv @@ -0,0 +1,122 @@ +PMTID p0 p1 p2 +332 17.9494 3.08921 -0.12284 +#334 6.65279 3.0573 -0.464513 +335 16.4242 3.08571 -0.124721 +336 13.721 3.08235 -0.13353 +#337 436.627 -31.3163 -36.6388 +338 9.01108 3.17324 -0.276284 +#339 6.13545 2.90173 -0.481356 +#340 6.07495 2.93731 -0.58905 +341 18.7845 3.09664 -0.131054 +343 12.7656 3.10667 -0.147391 +#344 6.12638 2.94238 -0.635735 +347 17.4949 3.10749 -0.149153 +348 17.8983 3.09234 -0.128647 +350 17.669 3.0891 -0.122669 +351 19.3944 3.08944 -0.120724 +353 8.42911 3.02959 0.102787 +354 15.5926 3.11266 -0.169246 +355 10.2884 3.09107 -0.185135 +356 7.64417 3.08218 0.168945 +357 16.2161 3.11699 -0.166343 +358 9.00471 3.08268 -0.171877 +359 14.4175 3.06029 -0.109384 +360 11.9095 3.10645 -0.184231 +361 8.17247 3.06593 -0.167536 +362 12.949 3.09306 -0.154836 +363 13.3929 3.09998 -0.158945 +364 10.4909 3.10289 -0.197985 +365 9.11282 3.0788 -0.185299 +366 7.43345 3.05244 0.15814 +367 14.3723 3.11048 -0.169921 +368 15.1638 3.12091 -0.179567 +369 7.85564 3.06208 -0.175745 +370 10.406 3.09483 -0.182604 +371 14.3153 3.11991 -0.184283 +372 20.8594 3.10104 -0.134912 +373 19.36 3.09268 -0.127957 +374 20.2707 3.10777 -0.144493 +375 19.384 3.0922 -0.125659 +376 17.9304 3.1021 -0.144849 +377 18.711 3.09601 -0.13753 +378 17.4338 3.1013 -0.144337 +379 20.6701 3.10214 -0.139123 +380 18.8799 3.08895 -0.123302 +381 18.0276 3.10058 -0.140683 +382 19.0526 3.16994 -0.212781 +383 19.2612 3.09715 -0.137064 +384 21.2208 3.10823 -0.14093 +385 20.2778 3.09937 -0.134291 +386 18.6713 3.09901 -0.137965 +387 19.5719 3.10468 -0.149018 +388 16.293 3.09538 -0.139371 +389 20.4549 3.10407 -0.137232 +390 17.1516 3.08933 -0.127157 +391 16.673 3.10151 -0.145288 +392 19.6285 3.10921 -0.145459 +393 24.4723 3.17837 -0.203509 +394 19.4412 3.10963 -0.147134 +395 22.5751 3.10263 -0.137148 +396 21.0331 3.10057 -0.138055 +397 19.1849 3.10283 -0.142099 +398 18.0316 3.09178 -0.129804 +399 18.707 3.10951 -0.150202 +400 18.1845 3.09103 -0.126774 +401 17.1005 3.10028 -0.140015 +402 17.3026 3.0921 -0.130513 +403 14.9342 3.08219 -0.122969 +404 17.604 3.17587 -0.216896 +405 11.9996 3.14612 -0.209651 +406 18.1137 3.11116 -0.154772 +407 18.8641 3.09702 -0.130414 +409 18.4361 3.09934 -0.140618 +410 17.2315 3.10588 -0.149666 +411 15.0428 3.10058 -0.145661 +412 19.0307 3.10149 -0.138611 +413 19.0157 3.09348 -0.12776 +414 20.8107 3.09968 -0.132717 +415 16.9133 3.08793 -0.12344 +417 19.8372 3.06649 -0.112973 +418 21.1259 3.0529 -0.11648 +419 18.1834 3.06772 -0.107166 +#420 6.77937 3.09534 -0.685827 +421 18.8051 3.06274 -0.115982 +422 19.4887 3.05639 -0.115288 +423 20.4322 3.05482 -0.111278 +424 20.5963 3.06444 -0.114702 +425 17.2319 3.07462 -0.121441 +426 28.5265 3.07966 -0.105101 +427 25.0327 3.08312 -0.105622 +428 17.9017 3.06745 -0.120512 +429 25.8197 3.07913 -0.112184 +430 21.8404 3.054 -0.115126 +432 18.5381 3.06253 -0.112199 +433 25.4477 3.07794 -0.109646 +434 18.1456 3.06517 -0.114085 +435 18.8795 3.06793 -0.113804 +436 28.6324 3.07859 -0.110805 +437 28.2059 3.07443 -0.105552 +438 20.8562 3.0513 -0.113908 +439 19.723 3.06511 -0.112248 +440 20.3786 3.05984 -0.119775 +441 18.9227 3.06643 -0.114837 +442 29.254 3.07889 -0.11029 +443 19.4492 3.0582 -0.109194 +446 20.4639 3.06413 -0.116817 +447 18.2738 3.06926 -0.111886 +448 26.1159 3.07479 -0.109723 +449 23.4118 3.07323 -0.10384 +450 11.6356 3.14276 -0.212469 +451 19.8326 3.06132 -0.108066 +452 19.2226 3.07337 -0.119869 +453 26.5169 3.07548 -0.107318 +454 27.3175 3.08095 -0.106168 +455 28.1968 3.08006 -0.113578 +456 36.6341 3.08324 -0.119822 +457 20.613 3.07381 -0.119642 +458 19.9613 3.05692 -0.119287 +459 19.8541 3.06219 -0.110582 +460 20.2442 3.06425 -0.11341 +461 19.451 3.07085 -0.109227 +462 18.9147 3.06294 -0.115266 +463 20.5881 3.05945 -0.110933 diff --git a/configfiles/PMTWaveformSim/PhaseIIADCHitFinderConfig b/configfiles/PMTWaveformSim/PhaseIIADCHitFinderConfig new file mode 100644 index 000000000..d05287eff --- /dev/null +++ b/configfiles/PMTWaveformSim/PhaseIIADCHitFinderConfig @@ -0,0 +1,11 @@ +verbosity 0 + +UseLEDWaveforms 0 + +PulseFindingApproach threshold +PulseWindowType dynamic +DefaultADCThreshold 7 +DefaultThresholdType relative + +EventBuilding 0 +MCWaveforms 1 diff --git a/configfiles/PMTWaveformSim/README.md b/configfiles/PMTWaveformSim/README.md new file mode 100644 index 000000000..78d3e2015 --- /dev/null +++ b/configfiles/PMTWaveformSim/README.md @@ -0,0 +1,25 @@ +# Configure files + +*********************** +#Description +********************** + +Configure files are simple text files for passing variables to the Tools. + +Text files are read by the Store class (src/Store) and automatically assigned to an internal map for the relevant Tool to use. + + +************************ +#Usage +************************ + +Any line starting with a "#" will be ignored by the Store, as will blank lines. + +Variables should be stored one per line as follows: + + +Name Value #Comments + + +Note: Only one value is permitted per name and they are stored in a string stream and template cast back to the type given. + diff --git a/configfiles/PMTWaveformSim/ToolChainConfig b/configfiles/PMTWaveformSim/ToolChainConfig new file mode 100644 index 000000000..5c69603a2 --- /dev/null +++ b/configfiles/PMTWaveformSim/ToolChainConfig @@ -0,0 +1,26 @@ +#ToolChain dynamic setup file + +##### Runtime Parameters ##### +verbose 0 ## Verbosity level of ToolChain +error_level 0 # 0= do not exit, 1= exit on unhandled errors only, 2= exit on unhandled errors and handled errors +attempt_recover 1 ## 1= will attempt to finalise if an execute fails +remote_port 24002 +IO_Threads 1 ## Number of threads for network traffic (~ 1/Gbps) + +###### Logging ##### +log_mode Interactive # Interactive=cout , Remote= remote logging system "serservice_name Remote_Logging" , Local = local file log; +log_local_path ./log +log_service LogStore + + +###### Service discovery ##### Ignore these settings for local analysis +service_publish_sec -1 +service_kick_sec -1 + +##### Tools To Add ##### +Tools_File configfiles/PMTWaveformSim/ToolsConfig ## list of tools to run and their config files + +##### Run Type ##### +Inline 2 ## number of Execute steps in program, -1 infinite loop that is ended by user +Interactive 0 ## set to 1 if you want to run the code interactively + diff --git a/configfiles/PMTWaveformSim/ToolsConfig b/configfiles/PMTWaveformSim/ToolsConfig new file mode 100644 index 000000000..42261ac88 --- /dev/null +++ b/configfiles/PMTWaveformSim/ToolsConfig @@ -0,0 +1,4 @@ +LoadGeometry LoadGeometry configfiles/PMTWaveformSim/LoadGeometryConfig +LoadWCSim LoadWCSim configfiles/PMTWaveformSim/LoadWCSimConfig +PMTWaveformSim PMTWaveformSim configfiles/PMTWaveformSim/PMTWaveformSimConfig +PhaseIIADCHitFinder PhaseIIADCHitFinder configfiles/PMTWaveformSim/PhaseIIADCHitFinderConfig From bcd0134d0ee8d17dc482e9da62cc80254bdd7d74 Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Tue, 3 Dec 2024 11:49:44 -0600 Subject: [PATCH 2/9] Big improvement. Fixes log norm functional form, better random throws for fit parameters, general cleanup --- UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 69 +++-- UserTools/PMTWaveformSim/PMTWaveformSim.h | 18 +- .../PMTWaveformSim/PMTWaveformSimConfig | 4 +- .../PMTfittingparametersWhitespace.csv | 241 +++++++++--------- 4 files changed, 185 insertions(+), 147 deletions(-) diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp index 9241bdfca..1dde9c641 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp @@ -115,9 +115,12 @@ bool PMTWaveformSim::Execute() uint16_t start_clocktick = (hit_t0 > fPrewindow)? hit_t0 - fPrewindow : 0; uint16_t end_clocktick = start_clocktick + fReadoutWindow; + // Randomly Sample the PMT parameters for each MCHit + SampleFitParameters(PMTID); + // loop over clock ticks for (uint16_t clocktick = start_clocktick; clocktick <= end_clocktick; clocktick += 1) { - uint16_t sample = CustomLogNormalPulse(hit_t0, clocktick, hit_charge, PMTID); + uint16_t sample = CustomLogNormalPulse(hit_t0, clocktick, hit_charge); // check if this hit time has been recorded // either set it or add to it @@ -149,6 +152,7 @@ bool PMTWaveformSim::Execute() }// end loop over PMTs // Publish the waveforms to the ANNIEEvent store + std::cout << "PMTWaveformSim: Saving RawADCDataMC with size: " << RawADCDataMC.size() << std::endl; m_data->Stores.at("ANNIEEvent")->Set("RawADCDataMC", RawADCDataMC); m_data->Stores.at("ANNIEEvent")->Set("CalibratedADCData", CalADCDataMC); @@ -181,26 +185,31 @@ bool PMTWaveformSim::LoadPMTParameters() } int pmtid; - double p0, p1, p2; + double p0, p1, p2, u00, u10, u11, u20, u21, u22; + std::string comma; std::string line; std::getline(infile, line); // Skipping the header line while (std::getline(infile, line)) { if (infile.fail()) { - logmessage = "PMTWaveformSim: Error reading from CSV file: "; + logmessage = "PMTWaveformSim: Error using CSV file: "; logmessage += fPMTParameterFile + "!"; Log(logmessage, v_error, verbosity); + return false; } // Skip any commented lines - if(line.find("#")!=std::string::npos) continue; + if(line.rfind("#",0)!=std::string::npos) continue; // Turn the line into a stringstream to extract the values std::stringstream ss(line); - ss >> pmtid >> p0 >> p1 >> p2; - - pmtParameters[pmtid] = std::make_tuple(p0, p1, p2); + ss >> pmtid >> comma >> p0 >> comma >> p1 >> comma >> p2 >> comma + >> u00 >> comma + >> u10 >> comma >> u11 >> comma + >> u20 >> comma >> u21 >> comma >> u22; + + fPMTParamMap[pmtid] = {p0, p1, p2, u00, u10, u11, u20, u21, u22}; logmessage = "PMTWaveformSim: Loaded parameters for PMTID " + std::to_string(pmtid) + ": "; logmessage += "p0 = " + std::to_string(p0); @@ -215,28 +224,48 @@ bool PMTWaveformSim::LoadPMTParameters() } //------------------------------------------------------------------------------ -uint16_t PMTWaveformSim::CustomLogNormalPulse(uint16_t hit_t0, uint16_t clocktick, double hit_charge, int PMTID) +bool PMTWaveformSim::SampleFitParameters(int pmtid) { - double p0, p1, p2; - if (pmtParameters.find(PMTID) != pmtParameters.end()) { - std::tie(p0, p1, p2) = pmtParameters[PMTID]; + PMTFitParams pmtParams; + if (fPMTParamMap.find(pmtid) != fPMTParamMap.end()) { + pmtParams = fPMTParamMap[pmtid]; } else { - logmessage = "PMTWaveformSim: PMTParameters not found for " + std::to_string(PMTID); + logmessage = "PMTWaveformSim: PMTParameters not found for " + std::to_string(pmtid); logmessage += ", using defaults: p0 = 17.49, p1 = 3.107, p2 = 0.1492"; Log(logmessage, v_warning, verbosity); - p0 = 17.49; - p1 = 3.107; - p2 = 0.1492; + // TODO make this a random sample as well + fP0 = 17.49; + fP1 = 3.107; + fP2 = 0.1492; + + return true; } + + // First sample a Gaussian with mean 0 and deviation 1 + double r0 = fRandom->Gaus(); + double r1 = fRandom->Gaus(); + double r2 = fRandom->Gaus(); + + // Convert to parameters that follow the fitted covariance matrix + fP0 = r0*pmtParams.u00 + pmtParams.p0; + fP1 = r0*pmtParams.u10 + r1*pmtParams.u11 + pmtParams.p1; + fP2 = r0*pmtParams.u20 + r1*pmtParams.u21 + r2*pmtParams.u22 + pmtParams.p2; - // The fit was performed in time units of ns, but we want to grab samples in clock ticks + return true; +} + +//------------------------------------------------------------------------------ +uint16_t PMTWaveformSim::CustomLogNormalPulse(uint16_t hit_t0, uint16_t clocktick, double hit_charge) +{ + //p0*exp( -0.5 * (log(x/p1)/p2)^2) + + // The fit was performed in time units of ns, but we pass samples in clock ticks double x = (double(clocktick) + fT0Offset - hit_t0) * NS_PER_ADC_SAMPLE; - // std::cout << " " << clocktick << ", " << x << std::endl; - double numerator = pow(log(x) - p1, 2); - double denom = (2 * pow(p2, 2)); - double amplitude = p0 * exp(-numerator / denom) * hit_charge; + double numerator = pow(log(x/fP1), 2); + double denom = (pow(fP2, 2)); + double amplitude = fP0 * exp(-0.5 * numerator/denom) * hit_charge; // Clip at 4095 and digitize to an integer return uint16_t((amplitude > 4095) ? 4095 : amplitude); diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.h b/UserTools/PMTWaveformSim/PMTWaveformSim.h index c46d8bf75..8008cadc1 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.h +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.h @@ -13,6 +13,14 @@ #include "TFile.h" #include "TRandom3.h" +struct PMTFitParams +{ + double p0; double p1; double p2; + double u00; + double u10; double u11; + double u20; double u21; double u22; +}; + /** * \class PMTWaveformSim @@ -35,7 +43,8 @@ class PMTWaveformSim: public Tool { bool LoadFromStores(); bool LoadPMTParameters(); - uint16_t CustomLogNormalPulse(uint16_t hit_t0, uint16_t t0_clocktick, double hit_charge, int PMTID); + bool SampleFitParameters(int pmtid); + uint16_t CustomLogNormalPulse(uint16_t hit_t0, uint16_t t0_clocktick, double hit_charge); void ConvertMapToWaveforms(const std::map &sample_map, const std::map> & parent_map, std::vector> &rawWaveforms, @@ -48,7 +57,7 @@ class PMTWaveformSim: public Tool { // To load from the ANNIEEvent std::map> *fMCHits = nullptr; - std::map> pmtParameters; + Geometry *fGeo = nullptr; // Config variables @@ -58,7 +67,10 @@ class PMTWaveformSim: public Tool { std::string fPMTParameterFile; TRandom3 *fRandom; - + + std::map fPMTParamMap; + double fP0, fP1, fP2; + bool fDebug; TFile *fOutFile; int fEvtNum = 0; diff --git a/configfiles/PMTWaveformSim/PMTWaveformSimConfig b/configfiles/PMTWaveformSim/PMTWaveformSimConfig index a271dfa82..d01603ac0 100644 --- a/configfiles/PMTWaveformSim/PMTWaveformSimConfig +++ b/configfiles/PMTWaveformSim/PMTWaveformSimConfig @@ -2,5 +2,5 @@ verbosity 0 PMTParameterFile configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv Prewindow 10 ReadoutWindow 35 -T0Offset 7 -MakeDebugFile 0 +T0Offset 0 +MakeDebugFile 1 diff --git a/configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv b/configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv index 5d42d4eda..76d5ec986 100644 --- a/configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv +++ b/configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv @@ -1,122 +1,119 @@ -PMTID p0 p1 p2 -332 17.9494 3.08921 -0.12284 -#334 6.65279 3.0573 -0.464513 -335 16.4242 3.08571 -0.124721 -336 13.721 3.08235 -0.13353 -#337 436.627 -31.3163 -36.6388 -338 9.01108 3.17324 -0.276284 -#339 6.13545 2.90173 -0.481356 -#340 6.07495 2.93731 -0.58905 -341 18.7845 3.09664 -0.131054 -343 12.7656 3.10667 -0.147391 -#344 6.12638 2.94238 -0.635735 -347 17.4949 3.10749 -0.149153 -348 17.8983 3.09234 -0.128647 -350 17.669 3.0891 -0.122669 -351 19.3944 3.08944 -0.120724 -353 8.42911 3.02959 0.102787 -354 15.5926 3.11266 -0.169246 -355 10.2884 3.09107 -0.185135 -356 7.64417 3.08218 0.168945 -357 16.2161 3.11699 -0.166343 -358 9.00471 3.08268 -0.171877 -359 14.4175 3.06029 -0.109384 -360 11.9095 3.10645 -0.184231 -361 8.17247 3.06593 -0.167536 -362 12.949 3.09306 -0.154836 -363 13.3929 3.09998 -0.158945 -364 10.4909 3.10289 -0.197985 -365 9.11282 3.0788 -0.185299 -366 7.43345 3.05244 0.15814 -367 14.3723 3.11048 -0.169921 -368 15.1638 3.12091 -0.179567 -369 7.85564 3.06208 -0.175745 -370 10.406 3.09483 -0.182604 -371 14.3153 3.11991 -0.184283 -372 20.8594 3.10104 -0.134912 -373 19.36 3.09268 -0.127957 -374 20.2707 3.10777 -0.144493 -375 19.384 3.0922 -0.125659 -376 17.9304 3.1021 -0.144849 -377 18.711 3.09601 -0.13753 -378 17.4338 3.1013 -0.144337 -379 20.6701 3.10214 -0.139123 -380 18.8799 3.08895 -0.123302 -381 18.0276 3.10058 -0.140683 -382 19.0526 3.16994 -0.212781 -383 19.2612 3.09715 -0.137064 -384 21.2208 3.10823 -0.14093 -385 20.2778 3.09937 -0.134291 -386 18.6713 3.09901 -0.137965 -387 19.5719 3.10468 -0.149018 -388 16.293 3.09538 -0.139371 -389 20.4549 3.10407 -0.137232 -390 17.1516 3.08933 -0.127157 -391 16.673 3.10151 -0.145288 -392 19.6285 3.10921 -0.145459 -393 24.4723 3.17837 -0.203509 -394 19.4412 3.10963 -0.147134 -395 22.5751 3.10263 -0.137148 -396 21.0331 3.10057 -0.138055 -397 19.1849 3.10283 -0.142099 -398 18.0316 3.09178 -0.129804 -399 18.707 3.10951 -0.150202 -400 18.1845 3.09103 -0.126774 -401 17.1005 3.10028 -0.140015 -402 17.3026 3.0921 -0.130513 -403 14.9342 3.08219 -0.122969 -404 17.604 3.17587 -0.216896 -405 11.9996 3.14612 -0.209651 -406 18.1137 3.11116 -0.154772 -407 18.8641 3.09702 -0.130414 -409 18.4361 3.09934 -0.140618 -410 17.2315 3.10588 -0.149666 -411 15.0428 3.10058 -0.145661 -412 19.0307 3.10149 -0.138611 -413 19.0157 3.09348 -0.12776 -414 20.8107 3.09968 -0.132717 -415 16.9133 3.08793 -0.12344 -417 19.8372 3.06649 -0.112973 -418 21.1259 3.0529 -0.11648 -419 18.1834 3.06772 -0.107166 -#420 6.77937 3.09534 -0.685827 -421 18.8051 3.06274 -0.115982 -422 19.4887 3.05639 -0.115288 -423 20.4322 3.05482 -0.111278 -424 20.5963 3.06444 -0.114702 -425 17.2319 3.07462 -0.121441 -426 28.5265 3.07966 -0.105101 -427 25.0327 3.08312 -0.105622 -428 17.9017 3.06745 -0.120512 -429 25.8197 3.07913 -0.112184 -430 21.8404 3.054 -0.115126 -432 18.5381 3.06253 -0.112199 -433 25.4477 3.07794 -0.109646 -434 18.1456 3.06517 -0.114085 -435 18.8795 3.06793 -0.113804 -436 28.6324 3.07859 -0.110805 -437 28.2059 3.07443 -0.105552 -438 20.8562 3.0513 -0.113908 -439 19.723 3.06511 -0.112248 -440 20.3786 3.05984 -0.119775 -441 18.9227 3.06643 -0.114837 -442 29.254 3.07889 -0.11029 -443 19.4492 3.0582 -0.109194 -446 20.4639 3.06413 -0.116817 -447 18.2738 3.06926 -0.111886 -448 26.1159 3.07479 -0.109723 -449 23.4118 3.07323 -0.10384 -450 11.6356 3.14276 -0.212469 -451 19.8326 3.06132 -0.108066 -452 19.2226 3.07337 -0.119869 -453 26.5169 3.07548 -0.107318 -454 27.3175 3.08095 -0.106168 -455 28.1968 3.08006 -0.113578 -456 36.6341 3.08324 -0.119822 -457 20.613 3.07381 -0.119642 -458 19.9613 3.05692 -0.119287 -459 19.8541 3.06219 -0.110582 -460 20.2442 3.06425 -0.11341 -461 19.451 3.07085 -0.109227 -462 18.9147 3.06294 -0.115266 -463 20.5881 3.05945 -0.110933 +PMT, p0, p1, p2, U00, U10, U11, U20, U21, U22, +332, 16.1815, 11.7835, 0.224469, 1.09912, 0.0012609, 0.20212, -0.0117725, -0.00387271, 0.0145493 +334, 7.92607, 10.2051, 0.657054, 0.528877, 0.00573739, 0.597992, -0.0455871, -0.0312974, 0.0506679 +335, 14.1205, 12.3273, 0.314469, 2.29598, -0.150105, 0.621752, -0.0458151, -0.00306655, 0.0480224 +336, 13.1547, 11.5986, 0.256954, 1.0054, 0.00256062, 0.250369, -0.0162861, -0.00559374, 0.0188073 +338, 8.63077, 12.4335, 0.474732, 0.942343, -0.131487, 0.755961, -0.0435303, -0.0117456, 0.0557825 +339, 6.5578, 8.94932, 0.585594, 0.803283, 0.109556, 0.91546, -0.0648363, -0.0612419, 0.068101 +340, 10.5316, 11.852, 0.358028, 0.937849, -0.0359459, 0.410756, -0.0258826, -0.00849393, 0.0305087 +341, 17.2763, 12.0659, 0.237495, 1.2176, -0.00123625, 0.222479, -0.0135174, -0.00415835, 0.0160533 +344, 11.145, 11.6382, 0.270989, 2.01325, -0.0741067, 0.521478, -0.0524951, -0.00598261, 0.0463826 +347, 17.066, 12.3058, 0.252206, 0.52976, -0.000135767, 0.111231, -0.00561469, -0.00219789, 0.00733403 +348, 16.6618, 11.9753, 0.235212, 1.21115, 0.0036329, 0.227827, -0.013784, -0.00467008, 0.0163823 +350, 17.4545, 11.778, 0.208261, 0.992566, 0.00912863, 0.161067, -0.00902673, -0.00359592, 0.0112417 +351, 18.4885, 11.7882, 0.196192, 1.21593, 0.0026219, 0.170627, -0.00999694, -0.00300509, 0.0123515 +353, 7.89973, 9.71231, 0.155955, 1.41864, -0.0200549, 0.332162, -0.0146225, -0.00466386, 0.0245054 +354, 14.7095, 12.3313, 0.298, 0.962958, -0.0163942, 0.26598, -0.0149568, -0.00475438, 0.0184653 +355, 10.255, 11.6652, 0.341621, 0.596657, -0.0118375, 0.262788, -0.0149194, -0.00620281, 0.0186959 +356, 7.51859, 10.2028, 0.314475, 0.819164, 0.0179374, 0.403051, -0.0278596, -0.0135264, 0.0325317 +357, 14.79, 12.6458, 0.321597, 0.905763, -0.0193769, 0.278233, -0.0150576, -0.00505714, 0.018792 +358, 9.09274, 11.0971, 0.314893, 0.706831, -0.0231053, 0.299926, -0.0182626, -0.00595521, 0.0228697 +360, 11.8003, 11.9891, 0.33628, 0.672946, -0.0155954, 0.258834, -0.0144044, -0.00548394, 0.0181108 +361, 8.3358, 10.5122, 0.329485, 0.582343, -0.00468082, 0.274878, -0.01737, -0.00772567, 0.0214473 +362, 12.5612, 11.7437, 0.291085, 0.815987, -0.00616299, 0.245508, -0.0149959, -0.00529733, 0.0179324 +363, 12.808, 11.9107, 0.289881, 0.692459, -0.0070296, 0.20929, -0.0118531, -0.00427004, 0.0147672 +364, 10.4694, 11.5381, 0.349791, 0.45393, 0.00594671, 0.209203, -0.0113758, -0.00645827, 0.0141728 +365, 8.96467, 10.8896, 0.38109, 0.509903, 0.00752423, 0.28113, -0.0167097, -0.00980058, 0.0202394 +366, 7.54604, 10.3277, 0.320294, 0.76851, -0.0255958, 0.358688, -0.0268422, -0.00812207, 0.030787 +367, 13.4516, 12.2728, 0.315725, 0.726559, -0.0103248, 0.235225, -0.0129952, -0.00486405, 0.0161507 +368, 13.7834, 12.4673, 0.323781, 0.625156, -0.00136558, 0.211405, -0.0112286, -0.00505615, 0.0139609 +369, 8.02919, 10.5823, 0.360006, 0.4664, -0.00974797, 0.257012, -0.0146779, -0.00726211, 0.0192091 +370, 10.2022, 11.356, 0.330295, 0.554959, 0.00188158, 0.236099, -0.0136587, -0.00666618, 0.0168183 +371, 13.6898, 12.3247, 0.324826, 0.638444, -0.00147777, 0.215787, -0.0115014, -0.00523692, 0.0143534 +372, 18.6865, 11.9436, 0.226671, 0.945701, -0.00285972, 0.152502, -0.00885944, -0.00257133, 0.010953 +373, 17.2248, 11.7201, 0.214925, 1.19781, 0.0196466, 0.208605, -0.0112031, -0.00547322, 0.0140196 +374, 17.653, 11.9499, 0.225296, 0.692605, 0.000785657, 0.120219, -0.00660593, -0.00228146, 0.00838178 +375, 17.3329, 11.626, 0.205215, 0.953283, 0.00644134, 0.149391, -0.00876697, -0.00316025, 0.0107825 +376, 16.3353, 11.7771, 0.234094, 1.10946, -0.01189, 0.203196, -0.0126647, -0.00287989, 0.0152296 +377, 16.4613, 11.7395, 0.228804, 1.33221, -0.0104825, 0.231954, -0.015608, -0.00349413, 0.0178923 +378, 15.5484, 11.894, 0.248814, 0.998359, -0.00496052, 0.207982, -0.0128779, -0.00379268, 0.0152893 +379, 18.4251, 11.8365, 0.216831, 1.4908, -0.00562787, 0.221649, -0.0150193, -0.00348075, 0.0170231 +380, 16.8922, 11.7073, 0.213033, 0.811405, 0.0104526, 0.142227, -0.00746291, -0.00346849, 0.00958732 +381, 15.8199, 11.8887, 0.247267, 0.886527, -0.00300425, 0.1825, -0.0109067, -0.00342333, 0.0132167 +382, 12.1736, 11.5981, 0.288064, 0.64403, 0.00322319, 0.203865, -0.0114588, -0.00514149, 0.01429 +383, 17.3771, 11.7654, 0.225753, 1.0952, -0.00326426, 0.183473, -0.011478, -0.00314874, 0.0136821 +384, 18.4667, 12.0561, 0.228487, 0.714233, 0.00195128, 0.122002, -0.00658039, -0.00242593, 0.00836149 +385, 16.6384, 11.8575, 0.223512, 1.42805, -0.0194552, 0.241996, -0.0158588, -0.00274704, 0.0184799 +386, 16.531, 11.9013, 0.231433, 0.915276, 0.00100232, 0.17039, -0.0100904, -0.00331765, 0.0122502 +387, 17.8784, 12.0554, 0.237911, 0.894593, -0.0030315, 0.160264, -0.00912539, -0.00281532, 0.0113281 +388, 15.2968, 11.6789, 0.240434, 1.38132, -0.0312206, 0.266723, -0.0179417, -0.0025785, 0.0207715 +389, 18.6947, 11.9855, 0.225579, 0.785306, 6.95106e-05, 0.128372, -0.00716647, -0.00236321, 0.0090084 +390, 16.193, 11.6874, 0.214387, 1.07159, -0.0034306, 0.182333, -0.0113206, -0.00296363, 0.0136509 +391, 14.3684, 11.5356, 0.239436, 0.885146, -0.00791972, 0.186319, -0.0116311, -0.00304685, 0.0140958 +392, 15.5213, 11.8151, 0.232703, 0.846511, -0.00272054, 0.166728, -0.00987177, -0.00296001, 0.0121191 +393, 14.4607, 11.9787, 0.276474, 0.851475, -0.0152815, 0.214403, -0.0123445, -0.0034262, 0.0153607 +394, 17.989, 12.0286, 0.230265, 0.800926, -0.00400728, 0.138462, -0.00763204, -0.00222846, 0.00972108 +395, 20.956, 11.916, 0.209918, 1.4117, -0.00115246, 0.184391, -0.0115401, -0.00307273, 0.0136481 +396, 19.4333, 12.0084, 0.2152, 1.20421, 0.00214268, 0.177116, -0.0107999, -0.00328542, 0.0128342 +397, 16.8615, 11.9109, 0.238467, 1.04534, -0.00730006, 0.192558, -0.011761, -0.00311293, 0.0141518 +398, 16.1561, 11.671, 0.224238, 0.987765, 0.00258986, 0.178236, -0.0109789, -0.00357651, 0.0131557 +399, 16.3963, 11.8413, 0.227827, 1.42694, -0.0291131, 0.24974, -0.0159717, -0.00210604, 0.0189558 +400, 17.1704, 11.7925, 0.214773, 1.01624, -0.00315399, 0.165391, -0.0100399, -0.00266275, 0.0122081 +401, 15.3736, 11.7643, 0.239775, 0.719101, -0.00404763, 0.147793, -0.00844082, -0.00256813, 0.0106377 +402, 16.4304, 11.5953, 0.213623, 0.891156, -0.0033035, 0.149204, -0.00896757, -0.00238849, 0.0110989 +403, 14.4771, 11.315, 0.214186, 0.770042, -3.11814e-05, 0.145043, -0.00872536, -0.00269165, 0.0108905 +404, 11.6303, 12.0387, 0.337009, 0.573103, -0.0126925, 0.226651, -0.0123472, -0.00486966, 0.0156723 +405, 7.82446, 10.8991, 0.339089, 0.707243, -0.0470227, 0.356936, -0.0235469, -0.00608166, 0.0285459 +406, 17.1667, 11.9454, 0.237854, 0.88653, -0.00666752, 0.162609, -0.00938048, -0.00256132, 0.0116844 +407, 17.7687, 11.9023, 0.212032, 0.766962, -0.000433398, 0.122873, -0.00689259, -0.00209822, 0.0087139 +409, 16.9692, 11.7387, 0.220713, 1.12227, -0.00308527, 0.187484, -0.0118391, -0.00317149, 0.0140611 +410, 14.9763, 11.9452, 0.264532, 0.964051, -0.0144364, 0.2222, -0.0132197, -0.00349243, 0.0161554 +411, 14.3011, 11.5125, 0.233957, 0.980853, -0.0170196, 0.199695, -0.0126145, -0.00239796, 0.0153161 +412, 17.4909, 11.8264, 0.216318, 1.03942, -0.00252854, 0.167459, -0.010243, -0.00276921, 0.012358 +413, 18.4095, 11.7611, 0.208316, 0.9201, 0.00638486, 0.140023, -0.00802475, -0.00299365, 0.00992059 +414, 19.0506, 11.7197, 0.200956, 1.04819, 0.00567565, 0.146363, -0.00868709, -0.00295905, 0.0105898 +415, 16.429, 11.5389, 0.207206, 0.924626, -0.00151963, 0.150925, -0.00887943, -0.0025144, 0.0111488 +417, 16.6145, 11.5931, 0.219003, 1.91434, -0.0116268, 0.303973, -0.0228075, -0.00461274, 0.0245922 +418, 18.1754, 11.5369, 0.1972, 2.26999, 0.00690851, 0.308702, -0.0209887, -0.00577008, 0.023937 +419, 16.7176, 11.4496, 0.197993, 1.79447, 0.0115074, 0.264508, -0.0184788, -0.0054933, 0.0206809 +420, 16.6491, 11.333, 0.185615, 5.22613, -0.0625622, 0.594426, -0.066662, -0.00452317, 0.0561069 +421, 16.5196, 11.5699, 0.206677, 1.66467, 0.00636307, 0.263751, -0.017533, -0.00517831, 0.020186 +422, 17.5206, 11.6951, 0.20264, 2.09836, 0.00128097, 0.300455, -0.0216512, -0.00519984, 0.0236865 +423, 17.3143, 11.7567, 0.205548, 1.74809, 0.00302894, 0.266245, -0.0173402, -0.00481607, 0.020097 +424, 18.8025, 11.6723, 0.192618, 2.25663, 0.00946812, 0.290165, -0.0202556, -0.00548813, 0.022506 +425, 15.4372, 11.6983, 0.246771, 2.1073, -0.0447784, 0.386148, -0.0326742, -0.00422519, 0.0324745 +426, 21.4214, 11.7092, 0.162415, 2.10618, 0.0157944, 0.216375, -0.0115835, -0.00436897, 0.0148912 +427, 18.3197, 11.5901, 0.185997, 1.94919, 0.0205332, 0.255456, -0.0166436, -0.00578503, 0.0190691 +428, 15.5427, 11.7414, 0.228688, 1.81489, -0.0129791, 0.322544, -0.0245963, -0.0050347, 0.0259873 +429, 19.8808, 11.5239, 0.169535, 2.25705, 0.0277746, 0.254673, -0.0150648, -0.00606648, 0.0181274 +430, 18.5201, 11.6785, 0.192243, 1.97639, 0.0124438, 0.271381, -0.0160084, -0.00547, 0.019617 +432, 16.2194, 11.5495, 0.21082, 1.79304, -0.00449995, 0.285874, -0.0203824, -0.00471571, 0.0227028 +433, 20.3546, 11.5811, 0.169168, 1.77143, 0.0203762, 0.197991, -0.0111087, -0.00461118, 0.0137774 +434, 16.8308, 11.6776, 0.204703, 1.90792, -0.00138537, 0.289198, -0.0201872, -0.0048371, 0.0225977 +435, 16.1769, 11.5845, 0.212605, 1.83411, 0.00378384, 0.296374, -0.0216042, -0.00565352, 0.0235066 +436, 20.6799, 11.6408, 0.167326, 2.08367, 0.0235311, 0.226964, -0.0128406, -0.00525132, 0.0157995 +437, 21.9742, 11.7023, 0.170203, 2.23103, 0.015895, 0.229293, -0.0133747, -0.00465486, 0.0163782 +438, 17.244, 11.5206, 0.194149, 1.94546, 0.0164538, 0.278271, -0.0186331, -0.00603257, 0.0212209 +439, 17.598, 11.5839, 0.20636, 1.92805, 0.00501941, 0.281054, -0.0198236, -0.00534138, 0.0220233 +440, 17.244, 11.8343, 0.22008, 2.05148, -0.00242051, 0.324049, -0.0238332, -0.00568393, 0.02552 +441, 16.7511, 11.6316, 0.212069, 1.98284, 0.00564093, 0.310735, -0.0224549, -0.00603168, 0.0244694 +442, 20.8555, 11.5429, 0.16751, 1.79699, 0.0185962, 0.193175, -0.0107737, -0.00436848, 0.0134829 +443, 16.8946, 11.3716, 0.191988, 1.61059, 0.0200682, 0.237307, -0.0148217, -0.00573007, 0.017558 +446, 17.0102, 11.8694, 0.220794, 1.88812, -0.00140432, 0.306529, -0.0219894, -0.00545122, 0.0238593 +447, 15.6919, 11.671, 0.22151, 1.93993, -0.0192759, 0.324627, -0.0256407, -0.00439161, 0.0266903 +448, 20.2622, 11.6803, 0.171468, 1.59172, 0.0152281, 0.185896, -0.00947263, -0.00403246, 0.0124752 +449, 21.0339, 11.5684, 0.165837, 1.96438, 0.0238181, 0.208036, -0.0117683, -0.00499108, 0.0144516 +450, 26.2343, 11.4588, 0.14365, 3.19853, 0.0102589, 0.233484, -0.0108197, -0.00385624, 0.0154123 +451, 17.1911, 11.5185, 0.194249, 1.88346, 0.0205039, 0.271199, -0.0182804, -0.00626395, 0.0206045 +452, 15.4282, 11.6637, 0.250749, 2.10438, -0.0259987, 0.39473, -0.0334018, -0.00610021, 0.0330829 +453, 20.0478, 11.6385, 0.175812, 1.86005, 0.0192341, 0.217135, -0.0127603, -0.00488615, 0.015464 +454, 19.5356, 11.5196, 0.161747, 1.91231, 0.0219558, 0.211209, -0.0117655, -0.00481501, 0.0146956 +455, 19.3469, 11.6801, 0.180487, 2.33217, 0.0243363, 0.28719, -0.0175967, -0.00644664, 0.0207875 +456, 24.8544, 11.7885, 0.158929, 2.13574, 0.010196, 0.189276, -0.00922597, -0.00348487, 0.0125871 +457, 17.9635, 11.9083, 0.215111, 2.04839, 0.00328706, 0.307679, -0.022277, -0.0057057, 0.0239149 +458, 16.5816, 11.775, 0.216282, 1.91451, 0.00494447, 0.311505, -0.0225166, -0.00601306, 0.0243323 +459, 19.0233, 11.6407, 0.190107, 1.84959, 0.00795076, 0.234735, -0.015644, -0.00443225, 0.0179219 +460, 18.1027, 11.6885, 0.20582, 1.95369, 0.00432308, 0.278251, -0.0194696, -0.00516732, 0.0216305 +461, 17.2489, 11.6514, 0.212362, 1.70302, 0.00220343, 0.2621, -0.0182581, -0.00486691, 0.0203933 +462, 16.3098, 11.5501, 0.212058, 2.27866, -0.00407155, 0.345676, -0.0289145, -0.0058746, 0.0290091 +463, 18.1769, 11.479, 0.192802, 1.82787, 0.0127039, 0.24785, -0.0159712, -0.00520418, 0.018721 From eb60c0da893ba57572857fb7d9de996e23c8cece Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Tue, 3 Dec 2024 12:57:54 -0500 Subject: [PATCH 3/9] Some fits are bad. Instead use some average parameters for each PMT type with p0 scaled based on the gain data. --- ...Whitespace.csv => PMTWaveformLognormFit.csv} | 17 ++++++++++------- configfiles/PMTWaveformSim/PMTWaveformSimConfig | 4 ++-- 2 files changed, 12 insertions(+), 9 deletions(-) rename configfiles/PMTWaveformSim/{PMTfittingparametersWhitespace.csv => PMTWaveformLognormFit.csv} (94%) diff --git a/configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv b/configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv similarity index 94% rename from configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv rename to configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv index 76d5ec986..8ba882b55 100644 --- a/configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv +++ b/configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv @@ -1,23 +1,25 @@ PMT, p0, p1, p2, U00, U10, U11, U20, U21, U22, 332, 16.1815, 11.7835, 0.224469, 1.09912, 0.0012609, 0.20212, -0.0117725, -0.00387271, 0.0145493 -334, 7.92607, 10.2051, 0.657054, 0.528877, 0.00573739, 0.597992, -0.0455871, -0.0312974, 0.0506679 -335, 14.1205, 12.3273, 0.314469, 2.29598, -0.150105, 0.621752, -0.0458151, -0.00306655, 0.0480224 +334, 16.9673, 11.8651, 0.2489, 2.9381, 0, 0.2207, 0, 0, 0.0472 +335, 16.4180, 11.8651, 0.2489, 2.9381, 0, 0.2207, 0, 0, 0.0472 336, 13.1547, 11.5986, 0.256954, 1.0054, 0.00256062, 0.250369, -0.0162861, -0.00559374, 0.0188073 -338, 8.63077, 12.4335, 0.474732, 0.942343, -0.131487, 0.755961, -0.0435303, -0.0117456, 0.0557825 -339, 6.5578, 8.94932, 0.585594, 0.803283, 0.109556, 0.91546, -0.0648363, -0.0612419, 0.068101 +338, 16.5613, 11.8651, 0.2489, 2.9381, 0, 0.2207, 0, 0, 0.0472 +339, 12.3941, 11.8651, 0.2489, 2.9381, 0, 0.2207, 0, 0, 0.0472 340, 10.5316, 11.852, 0.358028, 0.937849, -0.0359459, 0.410756, -0.0258826, -0.00849393, 0.0305087 341, 17.2763, 12.0659, 0.237495, 1.2176, -0.00123625, 0.222479, -0.0135174, -0.00415835, 0.0160533 +343, 3.1881, 11.8651, 0.2489, 2.9381, 0, 0.2207, 0, 0, 0.0472 344, 11.145, 11.6382, 0.270989, 2.01325, -0.0741067, 0.521478, -0.0524951, -0.00598261, 0.0463826 347, 17.066, 12.3058, 0.252206, 0.52976, -0.000135767, 0.111231, -0.00561469, -0.00219789, 0.00733403 348, 16.6618, 11.9753, 0.235212, 1.21115, 0.0036329, 0.227827, -0.013784, -0.00467008, 0.0163823 350, 17.4545, 11.778, 0.208261, 0.992566, 0.00912863, 0.161067, -0.00902673, -0.00359592, 0.0112417 351, 18.4885, 11.7882, 0.196192, 1.21593, 0.0026219, 0.170627, -0.00999694, -0.00300509, 0.0123515 -353, 7.89973, 9.71231, 0.155955, 1.41864, -0.0200549, 0.332162, -0.0146225, -0.00466386, 0.0245054 +353, 3.7712, 11.6034, 0.3268, 2.4673, 0, 0.7434, 0, 0, 0.0241 354, 14.7095, 12.3313, 0.298, 0.962958, -0.0163942, 0.26598, -0.0149568, -0.00475438, 0.0184653 355, 10.255, 11.6652, 0.341621, 0.596657, -0.0118375, 0.262788, -0.0149194, -0.00620281, 0.0186959 -356, 7.51859, 10.2028, 0.314475, 0.819164, 0.0179374, 0.403051, -0.0278596, -0.0135264, 0.0325317 +356, 4.7437, 11.6034, 0.3268, 2.4673, 0, 0.7434, 0, 0, 0.024 357, 14.79, 12.6458, 0.321597, 0.905763, -0.0193769, 0.278233, -0.0150576, -0.00505714, 0.018792 358, 9.09274, 11.0971, 0.314893, 0.706831, -0.0231053, 0.299926, -0.0182626, -0.00595521, 0.0228697 +359, 2.7887, 11.6034, 0.3268, 2.4673, 0, 0.7434, 0, 0, 0.0241 360, 11.8003, 11.9891, 0.33628, 0.672946, -0.0155954, 0.258834, -0.0144044, -0.00548394, 0.0181108 361, 8.3358, 10.5122, 0.329485, 0.582343, -0.00468082, 0.274878, -0.01737, -0.00772567, 0.0214473 362, 12.5612, 11.7437, 0.291085, 0.815987, -0.00616299, 0.245508, -0.0149959, -0.00529733, 0.0179324 @@ -66,6 +68,7 @@ PMT, p0, p1, p2, U00, U10, U11, U20, U21, U22, 405, 7.82446, 10.8991, 0.339089, 0.707243, -0.0470227, 0.356936, -0.0235469, -0.00608166, 0.0285459 406, 17.1667, 11.9454, 0.237854, 0.88653, -0.00666752, 0.162609, -0.00938048, -0.00256132, 0.0116844 407, 17.7687, 11.9023, 0.212032, 0.766962, -0.000433398, 0.122873, -0.00689259, -0.00209822, 0.0087139 +408, 16.4457, 11.7988, 0.2254, 1.4540, 0, 0.1640, 0, 0, 0.0136 409, 16.9692, 11.7387, 0.220713, 1.12227, -0.00308527, 0.187484, -0.0118391, -0.00317149, 0.0140611 410, 14.9763, 11.9452, 0.264532, 0.964051, -0.0144364, 0.2222, -0.0132197, -0.00349243, 0.0161554 411, 14.3011, 11.5125, 0.233957, 0.980853, -0.0170196, 0.199695, -0.0126145, -0.00239796, 0.0153161 @@ -76,7 +79,7 @@ PMT, p0, p1, p2, U00, U10, U11, U20, U21, U22, 417, 16.6145, 11.5931, 0.219003, 1.91434, -0.0116268, 0.303973, -0.0228075, -0.00461274, 0.0245922 418, 18.1754, 11.5369, 0.1972, 2.26999, 0.00690851, 0.308702, -0.0209887, -0.00577008, 0.023937 419, 16.7176, 11.4496, 0.197993, 1.79447, 0.0115074, 0.264508, -0.0184788, -0.0054933, 0.0206809 -420, 16.6491, 11.333, 0.185615, 5.22613, -0.0625622, 0.594426, -0.066662, -0.00452317, 0.0561069 +420, 22.7669, 11.6322, 0.1959, 2.3697, 0, 0.1129, 0, 0, 0.0240 421, 16.5196, 11.5699, 0.206677, 1.66467, 0.00636307, 0.263751, -0.017533, -0.00517831, 0.020186 422, 17.5206, 11.6951, 0.20264, 2.09836, 0.00128097, 0.300455, -0.0216512, -0.00519984, 0.0236865 423, 17.3143, 11.7567, 0.205548, 1.74809, 0.00302894, 0.266245, -0.0173402, -0.00481607, 0.020097 diff --git a/configfiles/PMTWaveformSim/PMTWaveformSimConfig b/configfiles/PMTWaveformSim/PMTWaveformSimConfig index d01603ac0..1cfd33811 100644 --- a/configfiles/PMTWaveformSim/PMTWaveformSimConfig +++ b/configfiles/PMTWaveformSim/PMTWaveformSimConfig @@ -1,6 +1,6 @@ verbosity 0 -PMTParameterFile configfiles/PMTWaveformSim/PMTfittingparametersWhitespace.csv +PMTParameterFile configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv Prewindow 10 ReadoutWindow 35 T0Offset 0 -MakeDebugFile 1 +MakeDebugFile 0 From f6a4350d9e37ec4ac70cd3d825de32a1f889d8ba Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Tue, 3 Dec 2024 13:22:06 -0500 Subject: [PATCH 4/9] Adding an informative README --- UserTools/PMTWaveformSim/README.md | 40 ++++++++++++++++-------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/UserTools/PMTWaveformSim/README.md b/UserTools/PMTWaveformSim/README.md index 4142083da..ba457de91 100644 --- a/UserTools/PMTWaveformSim/README.md +++ b/UserTools/PMTWaveformSim/README.md @@ -1,20 +1,22 @@ -# PMTWaveformSim - -PMTWaveformSim - -## Data - -Describe any data formats PMTWaveformSim creates, destroys, changes, or analyzes. E.G. - -**RawLAPPDData** `map>>` -* Takes this data from the `ANNIEEvent` store and finds the number of peaks +# PMTWaveformSim - -## Configuration - -Describe any configuration variables for PMTWaveformSim. - -``` -param1 value1 -param2 value2 -``` +PMTWaveformSim generates PMT waveforms based on the simulated PMT hits. For each MCHit a lognorm waveform is sampled based on fit parameters extracted from SPE waveforms in data. The fit parameters are randomly sampled in a way that conserved the covariance between parameters and captures the random behavior of a PMT's response. A full extended readout window (70 us) is sampled and waveforms from overlapping hits are added. A random baseline between 300 and 350 ADC is added. Random noise is applied, where the noise sigma is sampled from a gaussian with mean 1 and std dev of 0.25. The baseline and noise envelope values are hardcoded at the moment. Finally, any ADC counts above 4095 are clipped to mimick saturation. + +## Data + +**RawADCDataMC** `std::map> >` +* The raw simulated waveforms. The key is PMT channel number, then inner vector should always have size == 1 + +**CalibratedADCData** `std::map> >` +* The "calibrated" simulated waveforms. The baseline and sigma are the true, randomly-generated values. The key is PMT channel number, then inner vector should always have size == 1 + + +## Configuration + +``` +PMTParameterFile configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv # file containing the fit parameters and triangular Cholesky decomposed covariance matrix +Prewindow 10 # number of clock ticks before the MC hit time to begin sampling the fit function +ReadoutWindow 35 # number of clock ticks around the MC hit time over which waveforms are sampled +T0Offset 0 # A timing offset (in clock ticks) that can be used to align the pulse start time +MakeDebugFile 0 # Produce a root file containing all the simulated waveforms +``` From d9fab3fee3c7d603eaa2ffd4149dafee5885ccea Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Tue, 3 Dec 2024 12:22:50 -0600 Subject: [PATCH 5/9] Removing debug print statements --- UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp index 1dde9c641..a827a90ee 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp @@ -90,8 +90,6 @@ bool PMTWaveformSim::Execute() for (auto mcHitsIt : *fMCHits) { // Loop over the hit PMTs int PMTID = mcHitsIt.first; - std::cout << "PMT: " << PMTID << std::endl; - std::vector mcHits = mcHitsIt.second; // Generate waveform samples from the MC hits @@ -109,8 +107,6 @@ bool PMTWaveformSim::Execute() uint16_t hit_t0 = uint16_t(mcHit.GetTime() / NS_PER_ADC_SAMPLE); double hit_charge = mcHit.GetCharge(); - std::cout << mcHit.GetTime() << ", " << hit_charge << std::endl; - // Set the readout window but don't allow negative times uint16_t start_clocktick = (hit_t0 > fPrewindow)? hit_t0 - fPrewindow : 0; uint16_t end_clocktick = start_clocktick + fReadoutWindow; @@ -152,7 +148,6 @@ bool PMTWaveformSim::Execute() }// end loop over PMTs // Publish the waveforms to the ANNIEEvent store - std::cout << "PMTWaveformSim: Saving RawADCDataMC with size: " << RawADCDataMC.size() << std::endl; m_data->Stores.at("ANNIEEvent")->Set("RawADCDataMC", RawADCDataMC); m_data->Stores.at("ANNIEEvent")->Set("CalibratedADCData", CalADCDataMC); From a81abc1ef47147bcb1ff3155a9d6a2c8fa83911d Mon Sep 17 00:00:00 2001 From: Andrew Sutton <75585895+atcsutton@users.noreply.github.com> Date: Tue, 3 Dec 2024 13:59:42 -0500 Subject: [PATCH 6/9] Delete configfiles/PMTWaveformSim/DummyToolConfig --- configfiles/PMTWaveformSim/DummyToolConfig | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 configfiles/PMTWaveformSim/DummyToolConfig diff --git a/configfiles/PMTWaveformSim/DummyToolConfig b/configfiles/PMTWaveformSim/DummyToolConfig deleted file mode 100644 index 95cad88d2..000000000 --- a/configfiles/PMTWaveformSim/DummyToolConfig +++ /dev/null @@ -1,3 +0,0 @@ -# Dummy config file - -verbose 2 \ No newline at end of file From 4da2d4263816df86b3dfa92361c6e3f611798bac Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Tue, 3 Dec 2024 14:28:53 -0600 Subject: [PATCH 7/9] Remove elements of sample map so the find operation runs more quickly. Also change TGraph pointer to an object. --- UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 32 +++++++++++++-------- UserTools/PMTWaveformSim/PMTWaveformSim.h | 4 +-- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp index a827a90ee..86de01862 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp @@ -183,8 +183,6 @@ bool PMTWaveformSim::LoadPMTParameters() double p0, p1, p2, u00, u10, u11, u20, u21, u22; std::string comma; std::string line; - std::getline(infile, line); // Skipping the header line - while (std::getline(infile, line)) { if (infile.fail()) { logmessage = "PMTWaveformSim: Error using CSV file: "; @@ -194,8 +192,11 @@ bool PMTWaveformSim::LoadPMTParameters() return false; } + // Skip the header line + if (line.find("PMT") != std::string::npos) continue; + // Skip any commented lines - if(line.rfind("#",0)!=std::string::npos) continue; + if(line.rfind("#",0) != std::string::npos) continue; // Turn the line into a stringstream to extract the values std::stringstream ss(line); @@ -267,8 +268,8 @@ uint16_t PMTWaveformSim::CustomLogNormalPulse(uint16_t hit_t0, uint16_t clocktic } //------------------------------------------------------------------------------ -void PMTWaveformSim::ConvertMapToWaveforms(const std::map &sample_map, - const std::map> & parent_map, +void PMTWaveformSim::ConvertMapToWaveforms(std::map &sample_map, + std::map> & parent_map, std::vector> &rawWaveforms, std::vector> &calWaveforms, double noiseSigma, int baseline) @@ -283,18 +284,25 @@ void PMTWaveformSim::ConvertMapToWaveforms(const std::map &s int sample = std::round(noise + baseline); - // look for this tick in the sample map and add it - if (sample_map.find(tick) != sample_map.end()) + // look for this tick in the sample map and add it + // then remove it from the map to make the next find faster + if (sample_map.find(tick) != sample_map.end()) { sample += sample_map.at(tick); - + sample_map.erase(tick); + } rawSamples.push_back((sample > 4095) ? 4095 : sample); calSamples.push_back((rawSamples.back() - baseline) * ADC_TO_VOLT); + // Look for parent info associated with each tick + // need to transfer a set into a vector + // then erase the tick from the map std::vector innerParents; - if (parent_map.find(tick) != parent_map.end()) + if (parent_map.find(tick) != parent_map.end()) { std::copy(parent_map.at(tick).begin(), parent_map.at(tick).begin(), std::back_inserter(innerParents)); + parent_map.erase(tick); + } else innerParents.push_back(-5); @@ -349,14 +357,14 @@ void PMTWaveformSim::FillDebugGraphs(const std::map samples = waveform.Samples(); - TGraph* grTemp = new TGraph(); + TGraph grTemp = TGraph(); double sampleX = waveform.GetStartTime() / NS_PER_ADC_SAMPLE; for(auto sample : samples) { - grTemp->AddPoint(sampleX, sample); + grTemp.AddPoint(sampleX, sample); ++sampleX; } - grTemp->Write(grName.c_str()); + grTemp.Write(grName.c_str()); }// end loop over waveforms }// end loop over PMTs } diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.h b/UserTools/PMTWaveformSim/PMTWaveformSim.h index 8008cadc1..713a12b63 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.h +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.h @@ -45,8 +45,8 @@ class PMTWaveformSim: public Tool { bool LoadPMTParameters(); bool SampleFitParameters(int pmtid); uint16_t CustomLogNormalPulse(uint16_t hit_t0, uint16_t t0_clocktick, double hit_charge); - void ConvertMapToWaveforms(const std::map &sample_map, - const std::map> & parent_map, + void ConvertMapToWaveforms(std::map &sample_map, + std::map> & parent_map, std::vector> &rawWaveforms, std::vector> &calWaveforms, double noiseSigma, int baseline); From e7a70153d9f2148d903128ed682fcc3a1e34142a Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Fri, 6 Dec 2024 12:55:45 -0600 Subject: [PATCH 8/9] Fixing a few things here. First, removing the MCWaveforms. Don't really need them to interface with BackTracker. Second, adding the stop_time to the ADCPulse (have version control and default value for backward compatibility). Third, integrating BackTracker. Fourth, allowing for an upstream tool to signal to PhaseIIADCHitFinder, ClusterFinder, and BackTracker that a given Execute step will be skipped. This occurs if there are no MCHits to process or if no good waveforms are produced, and prevents red herring errors from being thrown. --- DataModel/ADCPulse.cpp | 4 +- DataModel/ADCPulse.h | 14 +- DataModel/Waveform.h | 54 ---- UserTools/BackTracker/BackTracker.cpp | 235 +++++++++++++++--- UserTools/BackTracker/BackTracker.h | 19 +- UserTools/ClusterFinder/ClusterFinder.cpp | 10 + UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 148 ++++++----- UserTools/PMTWaveformSim/PMTWaveformSim.h | 9 +- .../PhaseIIADCHitFinder.cpp | 68 +++-- .../PMTWaveformSim/PMTWaveformLognormFit.csv | 8 +- 10 files changed, 355 insertions(+), 214 deletions(-) diff --git a/DataModel/ADCPulse.cpp b/DataModel/ADCPulse.cpp index 84a3e3247..7b85a2c53 100755 --- a/DataModel/ADCPulse.cpp +++ b/DataModel/ADCPulse.cpp @@ -7,9 +7,9 @@ ADCPulse::ADCPulse(int TubeId, double start_time, double peak_time, double baseline, double sigma_baseline, unsigned long area, unsigned short raw_amplitude, double calibrated_amplitude, - double charge) : Hit(TubeId, start_time, charge), + double charge, double stop_time) : Hit(TubeId, start_time, charge), start_time_(start_time), peak_time_(peak_time), baseline_(baseline), sigma_baseline_(sigma_baseline), raw_area_(area), - raw_amplitude_(raw_amplitude), calibrated_amplitude_(calibrated_amplitude) + raw_amplitude_(raw_amplitude), calibrated_amplitude_(calibrated_amplitude), stop_time_(stop_time) { } diff --git a/DataModel/ADCPulse.h b/DataModel/ADCPulse.h index 891e840d0..4176fbdf9 100755 --- a/DataModel/ADCPulse.h +++ b/DataModel/ADCPulse.h @@ -23,7 +23,7 @@ class ADCPulse : public Hit { ADCPulse(int TubeId, double start_time, double peak_time, double baseline, double sigma_baseline, unsigned long raw_area, unsigned short raw_amplitude, double calibrated_amplitude, - double charge); + double charge, double stop_time = 0); // @brief Returns the start time (ns) of the pulse relative to the // start of its minibuffer @@ -33,6 +33,10 @@ class ADCPulse : public Hit { // start of its minibuffer inline double peak_time() const { return peak_time_; } + // @brief Returns the stop time (ns) of the pulse relative to the + // start of its minibuffer + inline double stop_time() const { return stop_time_; } + // @brief Returns the approximate baseline (ADC) used to calibrate the // pulse inline double baseline() const { return baseline_; } @@ -71,16 +75,24 @@ class ADCPulse : public Hit { ar & raw_area_; ar & raw_amplitude_; ar & calibrated_amplitude_; + if (version > 0) + ar & stop_time_; } protected: double start_time_; // ns since beginning of minibuffer double peak_time_; // ns since beginning of minibuffer + double stop_time_; // ns since beginning of minibuffer double baseline_; // mean (ADC) double sigma_baseline_; // standard deviation (ADC) unsigned long raw_area_; // (ADC * samples) unsigned short raw_amplitude_; // ADC double calibrated_amplitude_; // V + }; + +// Need to increment the class version since we added time as a new variable +// the version number ensures backward compatibility when serializing +BOOST_CLASS_VERSION(ADCPulse, 1) diff --git a/DataModel/Waveform.h b/DataModel/Waveform.h index f564a0d8b..28d1b71a3 100644 --- a/DataModel/Waveform.h +++ b/DataModel/Waveform.h @@ -55,58 +55,4 @@ class Waveform : public SerialisableObject{ } }; -template -class MCWaveform : public Waveform { - - friend class boost::serialization::access; - - public: - MCWaveform() : Waveform(), fParents(std::vector>{}) {} - MCWaveform(double tsin, std::vector samplesin, std::vector> theparents) : Waveform(tsin, samplesin), fParents(theparents) {} - virtual ~MCWaveform(){}; - - inline const std::vector>* GetParents() { return &fParents; } - inline const std::vector>& Parents() { return fParents; } - inline std::vector ParentsAtSample(int sample) { return fParents[sample]; } - - inline void SetParents(std::vector> parentsin){ fParents = parentsin; } - - inline Waveform GetBaseWaveform(){ return *this; } - - // Override the base print to add in parent info - bool Print() { - int verbose=0; - cout << "StartTime : " << this->fStartTime << endl; - cout << "NSamples : " << this->fSamples.size() << endl; - if(verbose){ - cout << "Samples (parents) : {"; - for(int samplei=0; samplei < this->fSamples.size(); samplei++){ - cout << this->fSamples.at(samplei); - - if (this->fParents.size() == this->fSamples.size()) { - cout << "("; - for (auto parent : this->ParentsAtSample(samplei)) { - cout << parent; - if ((samplei+1) != this->ParentsAtSample(samplei).size()) - cout << ", "; - else - cout << ")"; - } - } - - if ((samplei+1) != this ->fSamples.size()) cout << ", "; - } - cout<<"}"<> fParents; - -}; - #endif diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index ce80ada9b..bf72fc03d 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -27,6 +27,13 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ Log(logmessage, v_error, verbosity); } + bool gotMCWaveforms = m_variables.Get("MCWaveforms", fMCWaveforms); + if (!gotMCWaveforms) { + fMCWaveforms = false; + logmessage = "BackTracker::Initialize: \"MCWaveforms\" not set in the config, defaulting to false"; + Log(logmessage, v_error, verbosity); + } + // Set up the pointers we're going to save. No need to // delete them at Finalize, the store will handle it @@ -42,9 +49,11 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ //------------------------------------------------------------------------------ bool BackTracker::Execute() { - if (!LoadFromStores()) - return false; - + int load_status = LoadFromStores(); + if (load_status == 0) return false; + if (load_status == 2) return true; + + fClusterToBestParticleID ->clear(); fClusterToBestParticlePDG->clear(); fClusterEfficiency ->clear(); @@ -52,25 +61,79 @@ bool BackTracker::Execute() fClusterTotalCharge ->clear(); fParticleToTankTotalCharge.clear(); + SumParticleTankCharge(); - - // Loop over the clusters and do the things - for (std::pair>&& apair : *fClusterMapMC) { - int prtId = -5; - int prtPdg = -5; - double eff = -5; - double pur = -5; - double totalCharge = 0; - - MatchMCParticle(apair.second, prtId, prtPdg, eff, pur, totalCharge); - - fClusterToBestParticleID ->emplace(apair.first, prtId); - fClusterToBestParticlePDG->emplace(apair.first, prtPdg); - fClusterEfficiency ->emplace(apair.first, eff); - fClusterPurity ->emplace(apair.first, pur); - fClusterTotalCharge ->emplace(apair.first, totalCharge); - } + if (fMCWaveforms) { // using clusters of Hits made from simulated PMT pulses + + // Produce the map from channel ID to pulse time to MCHit index + //std::cout << "BT: Calling MapPulsesToParentIdxs" << std::endl; + bool gotPulseMap = MapPulsesToParentIdxs(); + if (!gotPulseMap) { + logmessage = "BackTracker: No good pulse map."; + Log(logmessage, v_error, verbosity); + return false; + } + + // Loop over the clusters + for (std::pair>&& apair : *fClusterMap) { + // Create a vector of MCHits associated with the vector of Hits + std::vector mcHits; + for (auto hit : apair.second) { + int channel_key = hit.GetTubeId(); + double hitTime = hit.GetTime(); + + // Catches if something goes wrong + // First make sure that we have the PMT in the outer map + // Then make sure we have the hit time in the inner map + if (fMapChannelToPulseTimeToMCHitIdx.find(channel_key) == fMapChannelToPulseTimeToMCHitIdx.end()) { + std::cout << "BackTracker: No hit on this PMT: " << channel_key << std::endl; + return false; + } + if (fMapChannelToPulseTimeToMCHitIdx.at(channel_key).find(hitTime) == fMapChannelToPulseTimeToMCHitIdx.at(channel_key).end()) { + std::cout << "BackTracker: No hit on this PMT: " << channel_key << " at this time: " << hitTime << std::endl; + return false; + } + + std::vector mcHitIdxVec = fMapChannelToPulseTimeToMCHitIdx[channel_key][hitTime]; + for (auto mcHitIdx : mcHitIdxVec) + mcHits.push_back((fMCHitsMap->at(channel_key)).at(mcHitIdx)); + }// end loop over cluster hits + + int prtId = -5; + int prtPdg = -5; + double eff = -5; + double pur = -5; + double totalCharge = 0; + + MatchMCParticle(mcHits, prtId, prtPdg, eff, pur, totalCharge); + + fClusterToBestParticleID ->emplace(apair.first, prtId); + fClusterToBestParticlePDG->emplace(apair.first, prtPdg); + fClusterEfficiency ->emplace(apair.first, eff); + fClusterPurity ->emplace(apair.first, pur); + fClusterTotalCharge ->emplace(apair.first, totalCharge); + + m_data->Stores.at("ANNIEEvent")->Set("MapChannelToPulseTimeToMCHitIdx", fMapChannelToPulseTimeToMCHitIdx ); + }// end loop over the cluster map + } else { // using clusters of MCHits + // Loop over the MC clusters and do the things + for (std::pair>&& apair : *fClusterMapMC) { + int prtId = -5; + int prtPdg = -5; + double eff = -5; + double pur = -5; + double totalCharge = 0; + + MatchMCParticle(apair.second, prtId, prtPdg, eff, pur, totalCharge); + + fClusterToBestParticleID ->emplace(apair.first, prtId); + fClusterToBestParticlePDG->emplace(apair.first, prtPdg); + fClusterEfficiency ->emplace(apair.first, eff); + fClusterPurity ->emplace(apair.first, pur); + fClusterTotalCharge ->emplace(apair.first, totalCharge); + }// end loop over the MC cluster map + }// end if/else fMCWaveforms m_data->Stores.at("ANNIEEvent")->Set("ClusterToBestParticleID", fClusterToBestParticleID ); m_data->Stores.at("ANNIEEvent")->Set("ClusterToBestParticlePDG", fClusterToBestParticlePDG); @@ -91,8 +154,8 @@ bool BackTracker::Finalise() //------------------------------------------------------------------------------ void BackTracker::SumParticleTankCharge() { - for (auto mcHitsIt : *fMCHitsMap) { - std::vector mcHits = mcHitsIt.second; + for (auto apair : *fMCHitsMap) { + std::vector mcHits = apair.second; for (uint mcHitIdx = 0; mcHitIdx < mcHits.size(); ++mcHitIdx) { // technically a MCHit could have multiple parents, but they don't appear to in practice @@ -101,8 +164,8 @@ void BackTracker::SumParticleTankCharge() if (parentIdxs.size() != 1) continue; int particleId = -5; - for (auto it : *fMCParticleIndexMap) { - if (it.second == parentIdxs[0]) particleId = it.first; + for (auto bpair : *fMCParticleIndexMap) { + if (bpair.second == parentIdxs[0]) particleId = bpair.first; } if (particleId == -5) continue; @@ -133,8 +196,8 @@ void BackTracker::MatchMCParticle(std::vector const &mchits, int &prtId, } int particleId = -5; - for (auto it : *fMCParticleIndexMap) { - if (it.second == parentIdxs[0]) particleId = it.first; + for (auto apair : *fMCParticleIndexMap) { + if (apair.second == parentIdxs[0]) particleId = apair.first; } if (particleId == -5) continue; @@ -175,39 +238,129 @@ void BackTracker::MatchMCParticle(std::vector const &mchits, int &prtId, } //------------------------------------------------------------------------------ -bool BackTracker::LoadFromStores() +bool BackTracker::MapPulsesToParentIdxs() { - // Grab the stuff we need from the stores - bool goodMCClusters = m_data->CStore.Get("ClusterMapMC", fClusterMapMC); - if (!goodMCClusters) { - std::cerr<<"BackTracker: no ClusterMapMC in the CStore!"<> > adcPulseMap; + bool goodADCPulses = m_data->Stores.at("ANNIEEvent")->Get("RecoADCHits", adcPulseMap); + if (!goodADCPulses) { + logmessage = "BackTracker: no RecoADCHits in the ANNIEEvent!"; + Log(logmessage, v_error, verbosity); return false; } + // Loop over the ADCPulses and find MCHits that fall within the start and stop times + // also record the pulse time to match the MCHits to the reco Hits + // Pulses are indexed by PMT id then stacked in a two-deep vector + + for (auto apair: adcPulseMap) { + int channel_key = apair.first; + + std::map> mapHitTimeToParents; + bool goodPulses = false; + for (auto pulseVec : apair.second) { + for (auto pulse : pulseVec) { + goodPulses = true; + double pulseTime = pulse.peak_time(); + + // Record the hit index if it occurred within the pulse window + // If there is only one hit the no need to check the times + std::vector mcHits = fMCHitsMap->at(channel_key); + if (mcHits.size() == 1) + mapHitTimeToParents[pulseTime].push_back(0); + else { + for (uint mcHitIdx = 0; mcHitIdx < mcHits.size(); ++mcHitIdx) { + double hitTime = mcHits[mcHitIdx].GetTime(); + // The hit finding has to contend with noise so allow for a 10 ns + // slew in the pulse start time (I know it seems large) + if ( hitTime + 10 >= pulse.start_time() && hitTime <= pulse.stop_time()) + mapHitTimeToParents[pulseTime].push_back(mcHitIdx); + }// end loop over MCHits + } + }// end loop over inner pulse vector + }// end loop over outer pulse vector + + if (mapHitTimeToParents.size() == 0 && goodPulses) { + logmessage = "BackTracker::MapPulsesToParentIdxs: No MCHits match with this pulse! PMT channel: "; + logmessage += std::to_string(channel_key); + Log(logmessage, v_error, verbosity); + return false; + } + + fMapChannelToPulseTimeToMCHitIdx.emplace(channel_key, std::move(mapHitTimeToParents)); + }// end loop over pulse map + + return true; +} + +//------------------------------------------------------------------------------ +int BackTracker::LoadFromStores() +{ + // Grab the stuff we need from the stores + bool goodAnnieEvent = m_data->Stores.count("ANNIEEvent"); if (!goodAnnieEvent) { - std::cerr<<"BackTracker: no ANNIEEvent store!"<Stores.at("ANNIEEvent")->Get("SkipExecute", skip); + if (goodSkipStatus && skip) { + logmessage = "BackTracker: An upstream tool told me to skip this event."; + Log(logmessage, v_warning, verbosity); + return 2; + } + bool goodMCHits = m_data->Stores.at("ANNIEEvent")->Get("MCHits", fMCHitsMap); if (!goodMCHits) { - std::cerr<<"BackTracker: no MCHits in the ANNIEEvent!"<CStore.Get("ClusterMap", fClusterMap); + if (!goodClusters) { + logmessage = "BackTracker: no ClusterMap in the CStore!"; + Log(logmessage, v_error, verbosity); + return 0; + } + + bool goodRecoHits = m_data->Stores.at("ANNIEEvent")->Get("Hits", fRecoHitsMap); + if (!goodRecoHits) { + logmessage = "BackTracker: no Hits in the ANNIEEvent!"; + Log(logmessage, v_error, verbosity); + return 0; + } + + } else { + bool goodMCClusters = m_data->CStore.Get("ClusterMapMC", fClusterMapMC); + if (!goodMCClusters) { + logmessage = "BackTracker: no ClusterMapMC in the CStore!"; + Log(logmessage, v_error, verbosity); + return 0; + } + }// end if/else fMCWaveforms bool goodMCParticles = m_data->Stores.at("ANNIEEvent")->Get("MCParticles", fMCParticles); if (!goodMCParticles) { - std::cerr<<"BackTracker: no MCParticles in the ANNIEEvent!"<Stores.at("ANNIEEvent")->Get("TrackId_to_MCParticleIndex", fMCParticleIndexMap); if (!goodMCParticleIndexMap) { - std::cerr<<"BackTracker: no TrackId_to_MCParticleIndex in the ANNIEEvent!"< const &mchits, int &prtId, int &prtPdg, double &eff, double &pur, double &totalCharge); ///< The meat and potatoes + + bool MapPulsesToParentIdxs(); private: // Things we need to pull out of the store std::map> *fMCHitsMap = nullptr; ///< All of the MCHits keyed by channel number - std::map> *fClusterMapMC = nullptr; ///< Clusters that we will be linking MCParticles to + std::map> *fRecoHitsMap = nullptr; ///< All of the reco Hits keyed by channel number + std::map> *fClusterMapMC = nullptr; ///< MCClusters that we will be linking MCParticles to + std::map> *fClusterMap = nullptr; ///< Clusters that we will be linking MCParticles to std::vector *fMCParticles = nullptr; ///< The true particles from the event std::map *fMCParticleIndexMap = nullptr; ///< Map between the particle Id and it's position in MCParticles vector // We'll calculate this map from MCHit parent particle to the total charge deposited throughout the tank // technically a MCHit could have multiple parents, but they don't appear to in practice // the key is particle Id and value is total tank charge - std::map fParticleToTankTotalCharge; + std::map fParticleToTankTotalCharge; + + // Used in the case of MC Waveforms + // Outer map: PMT channel ID to inner map + // Inner map: pulse time to vector of MCHit indexes that contribute + // The peak times of ADCPulses should match the time associated with reco Hits + std::map>> fMapChannelToPulseTimeToMCHitIdx; // We'll save out maps between the local cluster time and // the ID and PDG of the particle that contributed the most energy @@ -57,6 +68,8 @@ class BackTracker: public Tool { std::map *fClusterPurity = nullptr; std::map *fClusterTotalCharge = nullptr; + bool fMCWaveforms; + /// \brief verbosity levels: if 'verbosity' < this level, the message type will be logged. int verbosity; int v_error=0; diff --git a/UserTools/ClusterFinder/ClusterFinder.cpp b/UserTools/ClusterFinder/ClusterFinder.cpp index 5c002b119..8ce1966b8 100644 --- a/UserTools/ClusterFinder/ClusterFinder.cpp +++ b/UserTools/ClusterFinder/ClusterFinder.cpp @@ -164,6 +164,16 @@ bool ClusterFinder::Execute(){ //---------------------------------------------------------------------------- //---------------get the members of the ANNIEEvent---------------------------- //---------------------------------------------------------------------------- + + // An upstream tool may opt to skip this execution stage + // For example the PMTWaveformSim tool will skip events with no MCHits or if + // no waveforms are produced. + bool skip = false; + bool got_skip_status = m_data->Stores["ANNIEEvent"]->Get("SkipExecute", skip); + if (got_skip_status && skip) { + Log("ClusterFinder: An upstream tool told me to skip this event.",v_warning,verbose); + return true; + } m_data->Stores["ANNIEEvent"]->Get("EventNumber", evnum); //m_data->Stores["ANNIEEvent"]->Get("BeamStatus", BeamStatus); diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp index 86de01862..83ad4c5c7 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp @@ -79,12 +79,20 @@ bool PMTWaveformSim::Initialise(std::string configfile, DataModel &data) //------------------------------------------------------------------------------ bool PMTWaveformSim::Execute() { - - if (!LoadFromStores()) + int load_status = LoadFromStores(); + if (load_status == 0) { + ++fEvtNum; return false; + } + if (load_status == 2) { + m_data->Stores.at("ANNIEEvent")->Set("SkipExecute", true); + ++fEvtNum; + return true; + } + // The container for the data that we'll put into the ANNIEEvent - std::map> > RawADCDataMC; + std::map> > RawADCDataMC; std::map> > CalADCDataMC; for (auto mcHitsIt : *fMCHits) { // Loop over the hit PMTs @@ -96,13 +104,14 @@ bool PMTWaveformSim::Execute() // samples from hits that are close in time will be added together // key is hit time in clock ticks, value is amplitude std::map sample_map; - std::map> parent_map; // use a set so each parent is recorded only once - - for (const auto& mcHit : mcHits) { // Loop through each MCHit in the vector - // convert PMT hit time to clock ticks and "digitize" by converting to an int - // skip negative hit times, what does that even mean if we're not using the smeared digit time? - if (mcHit.GetTime() < 0) continue; + for (const auto& mcHit : mcHits) {// Loop through each MCHit in the vector + // skip negative hit times, what does that even mean if we're not using the smeared digit time? + // skip hit times past 70 us since that's our longest readout + if (mcHit.GetTime() < 0) continue; + if (mcHit.GetTime() > 70000) continue; + + // convert PMT hit time to clock ticks and "digitize" by converting to an int // MCHit time is in ns, but we're going to sample in clock ticks uint16_t hit_t0 = uint16_t(mcHit.GetTime() / NS_PER_ADC_SAMPLE); double hit_charge = mcHit.GetCharge(); @@ -113,48 +122,57 @@ bool PMTWaveformSim::Execute() // Randomly Sample the PMT parameters for each MCHit SampleFitParameters(PMTID); - + + // if (mcHits.size() < 5) + //std::cout << "PMTWaveformSim: " << fEvtNum << ", " << hit_t0 << ", " << start_clocktick << ", " << end_clocktick << std::endl; + // loop over clock ticks - for (uint16_t clocktick = start_clocktick; clocktick <= end_clocktick; clocktick += 1) { + for (uint16_t clocktick = start_clocktick; clocktick <= end_clocktick; clocktick += 1) { uint16_t sample = CustomLogNormalPulse(hit_t0, clocktick, hit_charge); // check if this hit time has been recorded // either set it or add to it - if (sample_map.find(clocktick) == sample_map.end()) { + if (sample_map.find(clocktick) == sample_map.end()) sample_map[clocktick] = sample; - parent_map[clocktick] = std::set(mcHit.GetParents()->begin(), mcHit.GetParents()->end()); - } else { - sample_map[clocktick] += sample; - - // Put the parents into the set - for (uint idx = 0; idx < mcHit.GetParents()->size(); ++idx) - parent_map[clocktick].insert(mcHit.GetParents()->at(idx)); - } + else + sample_map[clocktick] += sample; }// end loop over clock ticks }// end loop over mcHits - + // If there are no samples for this PMT then no need to do the rest + if (sample_map.empty()) continue; + + // Set the noise envelope and baseline for this PMT // The noise std dev appears to be normally distributed around 1 with sigma 0.25 double noiseSigma = fRandom->Gaus(1, 0.25); int basline = fRandom->Uniform(300, 350); // convert the sample map into a vector of Waveforms and put them into the container - std::vector> rawWaveforms; + std::vector> rawWaveforms; std::vector> calWaveforms; - ConvertMapToWaveforms(sample_map, parent_map, rawWaveforms, calWaveforms, noiseSigma, basline); + ConvertMapToWaveforms(sample_map, rawWaveforms, calWaveforms, noiseSigma, basline); RawADCDataMC.emplace(PMTID, rawWaveforms); CalADCDataMC.emplace(PMTID, calWaveforms); }// end loop over PMTs - // Publish the waveforms to the ANNIEEvent store - m_data->Stores.at("ANNIEEvent")->Set("RawADCDataMC", RawADCDataMC); - m_data->Stores.at("ANNIEEvent")->Set("CalibratedADCData", CalADCDataMC); + // Publish the waveforms to the ANNIEEvent store if we have them + if (RawADCDataMC.size()) { + m_data->Stores.at("ANNIEEvent")->Set("RawADCDataMC", RawADCDataMC); + m_data->Stores.at("ANNIEEvent")->Set("CalibratedADCData", CalADCDataMC); + } else { + logmessage = "PMTWaveformSim: No waveforms produced. Skipping!"; + Log(logmessage, v_warning, verbosity); + m_data->Stores.at("ANNIEEvent")->Set("SkipExecute", true); + return true; + } + if (fDebug) FillDebugGraphs(RawADCDataMC); ++fEvtNum; + m_data->Stores.at("ANNIEEvent")->Set("SkipExecute", false); return true; } @@ -268,75 +286,69 @@ uint16_t PMTWaveformSim::CustomLogNormalPulse(uint16_t hit_t0, uint16_t clocktic } //------------------------------------------------------------------------------ -void PMTWaveformSim::ConvertMapToWaveforms(std::map &sample_map, - std::map> & parent_map, - std::vector> &rawWaveforms, +void PMTWaveformSim::ConvertMapToWaveforms(const std::map &sample_map, + std::vector> &rawWaveforms, std::vector> &calWaveforms, double noiseSigma, int baseline) { - // All MC has extended readout - std::vector rawSamples; - std::vector calSamples; - std::vector> parents; - for (uint16_t tick = 0; tick < 34993; ++tick) { + // Clear the output waveforms just in case + rawWaveforms.clear(); + calWaveforms.clear(); + + // All MC has extended readout, but it's time consuming to draw 35000 noise samples + // Instead, only sample up to the maximum tick time. Then start with all baseline and + // add noise only to relvant ticks. The intermediate space between pulses will have no noise + + uint16_t maxTick = sample_map.rbegin()->first; + std::vector rawSamples(maxTick+1, baseline); + std::vector calSamples(maxTick+1, 0); + for (auto sample_pair : sample_map) { + uint16_t tick = sample_pair.first; + // Generate noise for each sample based on the std dev of the noise envelope double noise = fRandom->Gaus(0, noiseSigma); - int sample = std::round(noise + baseline); - - // look for this tick in the sample map and add it - // then remove it from the map to make the next find faster - if (sample_map.find(tick) != sample_map.end()) { - sample += sample_map.at(tick); - sample_map.erase(tick); - } - rawSamples.push_back((sample > 4095) ? 4095 : sample); - calSamples.push_back((rawSamples.back() - baseline) * ADC_TO_VOLT); - - // Look for parent info associated with each tick - // need to transfer a set into a vector - // then erase the tick from the map - std::vector innerParents; - if (parent_map.find(tick) != parent_map.end()) { - std::copy(parent_map.at(tick).begin(), parent_map.at(tick).begin(), - std::back_inserter(innerParents)); - parent_map.erase(tick); - } - else - innerParents.push_back(-5); - - parents.push_back(innerParents); + sample += sample_pair.second; + + rawSamples[tick] = (sample > 4095) ? 4095 : sample; + calSamples[tick] = (rawSamples[tick] - baseline) * ADC_TO_VOLT; } - // The start time in data is a timestamp. Don't have that for MC so just set to 0. - rawWaveforms.emplace_back(0, rawSamples, parents); + // The start time in data is a trigger timestamp. Don't have that for MC so just set to 0. + rawWaveforms.emplace_back(0, rawSamples); calWaveforms.emplace_back(0, calSamples, baseline, noiseSigma); } //------------------------------------------------------------------------------ -bool PMTWaveformSim::LoadFromStores() +int PMTWaveformSim::LoadFromStores() { bool goodAnnieEvent = m_data->Stores.count("ANNIEEvent"); if (!goodAnnieEvent) { logmessage = "PMTWaveformSim: no ANNIEEvent store!"; Log(logmessage, v_error, verbosity); - return false; + return 0; } bool goodMCHits = m_data->Stores.at("ANNIEEvent")->Get("MCHits", fMCHits); if (!goodMCHits) { - logmessage = "PMTWaveformSim: no MCHits in the ANNIEEvent!"; + logmessage = "PMTWaveformSim: no MCHits in the ANNIEEvent! "; Log(logmessage, v_error, verbosity); - return false; + return 0; } + + if (fMCHits->empty()) { + logmessage = "PMTWaveformSim: The MCHits map is empty! Skipping!"; + Log(logmessage, v_warning, verbosity); + return 2; + } + - return true; - + return 1; } //------------------------------------------------------------------------------ -void PMTWaveformSim::FillDebugGraphs(const std::map> > &RawADCDataMC) +void PMTWaveformSim::FillDebugGraphs(const std::map> > &RawADCDataMC) { for (auto itpair : RawADCDataMC) { std::string chanString = std::to_string(itpair.first); @@ -358,7 +370,7 @@ void PMTWaveformSim::FillDebugGraphs(const std::map samples = waveform.Samples(); TGraph grTemp = TGraph(); - double sampleX = waveform.GetStartTime() / NS_PER_ADC_SAMPLE; + double sampleX = waveform.GetStartTime(); for(auto sample : samples) { grTemp.AddPoint(sampleX, sample); ++sampleX; diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.h b/UserTools/PMTWaveformSim/PMTWaveformSim.h index 713a12b63..50b5b6aa8 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.h +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.h @@ -40,18 +40,17 @@ class PMTWaveformSim: public Tool { bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools. bool Execute(); ///< Execute function used to perform Tool purpose. bool Finalise(); ///< Finalise function used to clean up resources. - bool LoadFromStores(); + int LoadFromStores(); bool LoadPMTParameters(); bool SampleFitParameters(int pmtid); uint16_t CustomLogNormalPulse(uint16_t hit_t0, uint16_t t0_clocktick, double hit_charge); - void ConvertMapToWaveforms(std::map &sample_map, - std::map> & parent_map, - std::vector> &rawWaveforms, + void ConvertMapToWaveforms(const std::map &sample_map, + std::vector> &rawWaveforms, std::vector> &calWaveforms, double noiseSigma, int baseline); - void FillDebugGraphs(const std::map> > &RawADCDataMC); + void FillDebugGraphs(const std::map> > &RawADCDataMC); private: diff --git a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp index cb05bd7d2..b6bc414b7 100644 --- a/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp +++ b/UserTools/PhaseIIADCHitFinder/PhaseIIADCHitFinder.cpp @@ -78,7 +78,9 @@ bool PhaseIIADCHitFinder::Initialise(std::string config_filename, DataModel& dat m_data->CStore.Get("AuxChannelNumToTypeMap",AuxChannelNumToTypeMap); // Get the timing offsets - m_data->CStore.Get("ChannelNumToTankPMTTimingOffsetMap",ChannelKeyToTimingOffsetMap); + // don't need them for MC + if (!mc_waveforms) + m_data->CStore.Get("ChannelNumToTankPMTTimingOffsetMap",ChannelKeyToTimingOffsetMap); //Recreate maps that were deleted with ANNIEEvent->Delete() ANNIEEventBuilder tool hit_map = new std::map>; @@ -137,30 +139,21 @@ bool PhaseIIADCHitFinder::Execute() { got_rawaux_data = annie_event->Get("RawADCAuxData", raw_aux_waveform_map); if (mc_waveforms) { - raw_waveform_map.clear(); - - std::map > > mc_waveform_map; - bool got_raw_mc = annie_event->Get("RawADCDataMC", mc_waveform_map); - if (!got_raw_mc) { - Log("Error: The PhaseIIADCHitFinder tool could not find the RawADCDataMC entry", v_error, - verbosity); - return false; - } + got_raw_data = annie_event->Get("RawADCDataMC", raw_waveform_map); + + // Some executes are skipped if there are no MCHits or waveforms produced + // Put the cleared maps into the ANNIEEvent to ensure that a downstream + // tool doesn't grab a map from a previous event + bool skip = false; + bool got_skip_status = annie_event->Get("SkipExecute", skip); + if (got_skip_status && skip) { + Log("PhaseIIADCHitFinder: An upstream tool told me to skip this event.",v_warning,verbosity); + + m_data->Stores.at("ANNIEEvent")->Set("RecoADCHits", pulse_map); + m_data->Stores.at("ANNIEEvent")->Set("RecoADCAuxHits", aux_pulse_map); - // Slice off the MC info to get a map of Waveforms - // Need to grab the parents for use in BackTracker later - for (auto itpair : mc_waveform_map) { - for (auto mc_waveform : itpair.second) - raw_waveform_map[itpair.first].push_back(mc_waveform.GetBaseWaveform()); + return true; } - - if (raw_waveform_map.size() == mc_waveform_map.size()) - got_raw_data = true; - else { - Log("Error: The PhaseIIADCHitFinder tool could not extract Waveforms from the MCWaveforms.", v_error, - verbosity); - return false; - } }// end if mc_waveforms } @@ -176,9 +169,9 @@ bool PhaseIIADCHitFinder::Execute() { return false; } else if ( raw_waveform_map.empty() ) { - Log("Error: The PhaseIIADCHitFinder tool found an empty RawADCData entry. Skipping.", v_error, + Log("Error: The PhaseIIADCHitFinder tool found an empty RawADCData entry.", v_error, verbosity); - return true; + return false; } // Load the maps containing the ADC calibrated waveform data @@ -205,9 +198,9 @@ bool PhaseIIADCHitFinder::Execute() { return false; } else if ( calibrated_waveform_map.empty() ) { - Log("Error: The PhaseIIADCHitFinder tool found an empty CalibratedADCData entry. Skipping", + Log("Error: The PhaseIIADCHitFinder tool found an empty CalibratedADCData entry.", v_error, verbosity); - return true; + return false; } //Find pulses in the raw detector data @@ -746,7 +739,7 @@ std::vector PhaseIIADCHitFinder::find_pulses_bywindow( if(it != ChannelKeyToTimingOffsetMap.end()){ //Timing offset is available timing_offset = ChannelKeyToTimingOffsetMap.at(channel_key); } else { - if(verbosity>2){ + if(verbosity>2 && !mc_waveforms){ std::cout << "Didn't find Timing offset for channel " << channel_key << std::endl; } } @@ -754,10 +747,11 @@ std::vector PhaseIIADCHitFinder::find_pulses_bywindow( // Store the freshly made pulse in the vector of found pulses pulses.emplace_back(channel_key, ( wmin * NS_PER_SAMPLE )-timing_offset, - (peak_sample * NS_PER_SAMPLE)-timing_offset, + ( peak_sample * NS_PER_SAMPLE )-timing_offset, calibrated_minibuffer_data.GetBaseline(), calibrated_minibuffer_data.GetSigmaBaseline(), - raw_area, max_ADC, calibrated_amplitude, charge); + raw_area, max_ADC, calibrated_amplitude, charge, + ( wmax * NS_PER_SAMPLE )-timing_offset); } return pulses; } @@ -860,7 +854,7 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( if(it != ChannelKeyToTimingOffsetMap.end()){ //Timing offset is available timing_offset = ChannelKeyToTimingOffsetMap.at(channel_key); } else { - if(verbosity>v_error){ + if(verbosity>v_error && !mc_waveforms){ std::cout << "PhaseIIADCHitFinder: Didn't find Timing offset for channel... setting this channel's offset to 0ns" << channel_key << std::endl; } } @@ -868,10 +862,11 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( // Store the freshly made pulse in the vector of found pulses pulses.emplace_back(channel_key, ( pulse_start_sample * NS_PER_SAMPLE )-timing_offset, - (peak_sample * NS_PER_SAMPLE)-timing_offset, + ( peak_sample * NS_PER_SAMPLE )-timing_offset, calibrated_minibuffer_data.GetBaseline(), calibrated_minibuffer_data.GetSigmaBaseline(), - raw_area, max_ADC, calibrated_amplitude, charge); + raw_area, max_ADC, calibrated_amplitude, charge, + ( pulse_end_sample * NS_PER_SAMPLE )-timing_offset); } @@ -937,7 +932,7 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( if(it != ChannelKeyToTimingOffsetMap.end()){ //Timing offset is available timing_offset = ChannelKeyToTimingOffsetMap.at(channel_key); } else { - if(verbosity>v_error){ + if(verbosity>v_error && !mc_waveforms){ std::cout << "PhaseIIADCHitFinder: Didn't find Timing offset for channel... setting this channel's offset to 0ns" << channel_key << std::endl; } } @@ -987,10 +982,11 @@ std::vector PhaseIIADCHitFinder::find_pulses_bythreshold( // Store the freshly made pulse in the vector of found pulses pulses.emplace_back(channel_key, ( pulse_start_sample * NS_PER_ADC_SAMPLE )-timing_offset, - (hit_time * NS_PER_ADC_SAMPLE)-timing_offset, // interpolated hit time + ( hit_time * NS_PER_ADC_SAMPLE )-timing_offset, // interpolated hit time calibrated_minibuffer_data.GetBaseline(), calibrated_minibuffer_data.GetSigmaBaseline(), - raw_area, max_ADC, calibrated_amplitude, charge); + raw_area, max_ADC, calibrated_amplitude, charge, + ( pulse_end_sample * NS_PER_ADC_SAMPLE )-timing_offset); } } } else { diff --git a/configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv b/configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv index 8ba882b55..0e7dc0969 100644 --- a/configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv +++ b/configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv @@ -7,19 +7,19 @@ PMT, p0, p1, p2, U00, U10, U11, U20, U21, U22, 339, 12.3941, 11.8651, 0.2489, 2.9381, 0, 0.2207, 0, 0, 0.0472 340, 10.5316, 11.852, 0.358028, 0.937849, -0.0359459, 0.410756, -0.0258826, -0.00849393, 0.0305087 341, 17.2763, 12.0659, 0.237495, 1.2176, -0.00123625, 0.222479, -0.0135174, -0.00415835, 0.0160533 -343, 3.1881, 11.8651, 0.2489, 2.9381, 0, 0.2207, 0, 0, 0.0472 +#343, 3.1881, 11.8651, 0.2489, 2.9381, 0, 0.2207, 0, 0, 0.0472 344, 11.145, 11.6382, 0.270989, 2.01325, -0.0741067, 0.521478, -0.0524951, -0.00598261, 0.0463826 347, 17.066, 12.3058, 0.252206, 0.52976, -0.000135767, 0.111231, -0.00561469, -0.00219789, 0.00733403 348, 16.6618, 11.9753, 0.235212, 1.21115, 0.0036329, 0.227827, -0.013784, -0.00467008, 0.0163823 350, 17.4545, 11.778, 0.208261, 0.992566, 0.00912863, 0.161067, -0.00902673, -0.00359592, 0.0112417 351, 18.4885, 11.7882, 0.196192, 1.21593, 0.0026219, 0.170627, -0.00999694, -0.00300509, 0.0123515 -353, 3.7712, 11.6034, 0.3268, 2.4673, 0, 0.7434, 0, 0, 0.0241 +#353, 3.7712, 11.6034, 0.3268, 2.4673, 0, 0.7434, 0, 0, 0.0241 354, 14.7095, 12.3313, 0.298, 0.962958, -0.0163942, 0.26598, -0.0149568, -0.00475438, 0.0184653 355, 10.255, 11.6652, 0.341621, 0.596657, -0.0118375, 0.262788, -0.0149194, -0.00620281, 0.0186959 -356, 4.7437, 11.6034, 0.3268, 2.4673, 0, 0.7434, 0, 0, 0.024 +#356, 4.7437, 11.6034, 0.3268, 2.4673, 0, 0.7434, 0, 0, 0.024 357, 14.79, 12.6458, 0.321597, 0.905763, -0.0193769, 0.278233, -0.0150576, -0.00505714, 0.018792 358, 9.09274, 11.0971, 0.314893, 0.706831, -0.0231053, 0.299926, -0.0182626, -0.00595521, 0.0228697 -359, 2.7887, 11.6034, 0.3268, 2.4673, 0, 0.7434, 0, 0, 0.0241 +#359, 2.7887, 11.6034, 0.3268, 2.4673, 0, 0.7434, 0, 0, 0.0241 360, 11.8003, 11.9891, 0.33628, 0.672946, -0.0155954, 0.258834, -0.0144044, -0.00548394, 0.0181108 361, 8.3358, 10.5122, 0.329485, 0.582343, -0.00468082, 0.274878, -0.01737, -0.00772567, 0.0214473 362, 12.5612, 11.7437, 0.291085, 0.815987, -0.00616299, 0.245508, -0.0149959, -0.00529733, 0.0179324 From d919fa7bfcde612eeed8c0fc05a57de342558ed0 Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Fri, 6 Dec 2024 13:43:47 -0600 Subject: [PATCH 9/9] Steven noticed a bug. I forgot to divide by the number of baslined when averaging them. --- UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp b/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp index 82c5d45fc..34b0fda62 100644 --- a/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp +++ b/UserTools/PhaseIIADCCalibrator/PhaseIIADCCalibrator.cpp @@ -780,6 +780,7 @@ PhaseIIADCCalibrator::make_calibrated_waveforms_ze3ra_multi( bl_estimates_mean += baselines[idx]; bl_estimates_sigma += pow(sigma_baselines[idx],2); } + bl_estimates_mean = bl_estimates_mean / baselines.size(); bl_estimates_sigma = sqrt(bl_estimates_sigma / double(baselines.size())); // ComputeMeanAndVariance(baselines, bl_estimates_mean, bl_estimates_var, std::numeric_limits::max(), 0, 7);