From eb71091a7f21019038959a00a4e6abc5596b472d Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Thu, 21 Nov 2024 13:12:56 -0600 Subject: [PATCH 01/15] 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 02/15] 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 03/15] 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 04/15] 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 05/15] 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 06/15] 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 07/15] 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 08/15] 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 09/15] 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); From 86dfa91e79ae314ab64c860d46b63605a3ccf070 Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Wed, 22 Jan 2025 09:14:57 -0600 Subject: [PATCH 10/15] Fix to the timing digitization. Now the raw hit time is passed to the pulse function. Also allowing the T0Offset parameter to be positive or negative. --- UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 29 +++++++++++---------- UserTools/PMTWaveformSim/PMTWaveformSim.h | 4 +-- UserTools/PMTWaveformSim/README.md | 2 +- 3 files changed, 18 insertions(+), 17 deletions(-) diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp index 83ad4c5c7..ceafffa5d 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp @@ -110,22 +110,19 @@ bool PMTWaveformSim::Execute() // 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); + + // Grab the hit time (also converted to clock ticks) and the charge + double hit_t0 = mcHit.GetTime(); + uint16_t t0_ticks = uint16_t(hit_t0 / NS_PER_ADC_SAMPLE); double hit_charge = mcHit.GetCharge(); - // Set the readout window but don't allow negative times - uint16_t start_clocktick = (hit_t0 > fPrewindow)? hit_t0 - fPrewindow : 0; + // Set the readout window in clock ticks, but don't allow negative times + uint16_t start_clocktick = (t0_ticks/ NS > fPrewindow)? t0_ticks - fPrewindow : 0; uint16_t end_clocktick = start_clocktick + fReadoutWindow; // 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) { uint16_t sample = CustomLogNormalPulse(hit_t0, clocktick, hit_charge); @@ -198,6 +195,9 @@ bool PMTWaveformSim::LoadPMTParameters() } int pmtid; + // Stored fit parameters. + // p0, p1, and p2 are the mean values of the lognorm + // the u's are for sampling in a way that follows the fit covariance double p0, p1, p2, u00, u10, u11, u20, u21, u22; std::string comma; std::string line; @@ -270,15 +270,16 @@ bool PMTWaveformSim::SampleFitParameters(int pmtid) } //------------------------------------------------------------------------------ -uint16_t PMTWaveformSim::CustomLogNormalPulse(uint16_t hit_t0, uint16_t clocktick, double hit_charge) +uint16_t PMTWaveformSim::CustomLogNormalPulse(double 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; - - double numerator = pow(log(x/fP1), 2); - double denom = (pow(fP2, 2)); + double x = (clocktick + fT0Offset) * NS_PER_ADC_SAMPLE - hit_t0; + + double numerator = log(x/fP1); + numerator = val * val; + double denom = fP2 * fP2; double amplitude = fP0 * exp(-0.5 * numerator/denom) * hit_charge; // Clip at 4095 and digitize to an integer diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.h b/UserTools/PMTWaveformSim/PMTWaveformSim.h index 50b5b6aa8..5fabf2ff8 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.h +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.h @@ -44,7 +44,7 @@ class PMTWaveformSim: public Tool { bool LoadPMTParameters(); bool SampleFitParameters(int pmtid); - uint16_t CustomLogNormalPulse(uint16_t hit_t0, uint16_t t0_clocktick, double hit_charge); + uint16_t CustomLogNormalPulse(double hit_t0, uint16_t t0_clocktick, double hit_charge); void ConvertMapToWaveforms(const std::map &sample_map, std::vector> &rawWaveforms, std::vector> &calWaveforms, @@ -62,7 +62,7 @@ class PMTWaveformSim: public Tool { // Config variables uint16_t fPrewindow; uint16_t fReadoutWindow; - uint16_t fT0Offset; + int fT0Offset; std::string fPMTParameterFile; TRandom3 *fRandom; diff --git a/UserTools/PMTWaveformSim/README.md b/UserTools/PMTWaveformSim/README.md index ba457de91..e4d4709a1 100644 --- a/UserTools/PMTWaveformSim/README.md +++ b/UserTools/PMTWaveformSim/README.md @@ -17,6 +17,6 @@ PMTWaveformSim generates PMT waveforms based on the simulated PMT hits. For each 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 +T0Offset 0 # A timing offset (in clock ticks) that can be used to align the pulse start time (can be positive or negative) MakeDebugFile 0 # Produce a root file containing all the simulated waveforms ``` From 3d94f7151aa92cd54ea4306ace594f71256c60a6 Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Wed, 22 Jan 2025 09:38:05 -0600 Subject: [PATCH 11/15] Fixing the event numbering for the debug plots so that they pull from the data store. --- UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 14 +++----------- UserTools/PMTWaveformSim/PMTWaveformSim.h | 1 - 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp index ceafffa5d..9e7293ac1 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp @@ -80,15 +80,6 @@ bool PMTWaveformSim::Initialise(std::string configfile, DataModel &data) bool PMTWaveformSim::Execute() { 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 @@ -168,7 +159,6 @@ bool PMTWaveformSim::Execute() if (fDebug) FillDebugGraphs(RawADCDataMC); - ++fEvtNum; m_data->Stores.at("ANNIEEvent")->Set("SkipExecute", false); return true; } @@ -366,7 +356,9 @@ void PMTWaveformSim::FillDebugGraphs(const std::mapStores.at("ANNIEEvent")->Get("EventNumber",evtNum); + std::string grName = ("wf_" + std::to_string(evtNum) + "_" + std::to_string(wfIdx)); // Make the graph std::vector samples = waveform.Samples(); diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.h b/UserTools/PMTWaveformSim/PMTWaveformSim.h index 5fabf2ff8..3158f8aa7 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.h +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.h @@ -72,7 +72,6 @@ class PMTWaveformSim: public Tool { bool fDebug; TFile *fOutFile; - int fEvtNum = 0; int verbosity; int v_error=0; From 69f76566e09d97887b749e84fafdf8705e73e9e0 Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Wed, 22 Jan 2025 10:01:34 -0600 Subject: [PATCH 12/15] Fixing dumb bugs. I got lazy before and didn't do a test build of my changes... --- UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp index 9e7293ac1..ed4466659 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.cpp +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.cpp @@ -108,7 +108,7 @@ bool PMTWaveformSim::Execute() double hit_charge = mcHit.GetCharge(); // Set the readout window in clock ticks, but don't allow negative times - uint16_t start_clocktick = (t0_ticks/ NS > fPrewindow)? t0_ticks - fPrewindow : 0; + uint16_t start_clocktick = (t0_ticks > fPrewindow)? t0_ticks - fPrewindow : 0; uint16_t end_clocktick = start_clocktick + fReadoutWindow; // Randomly Sample the PMT parameters for each MCHit @@ -268,7 +268,7 @@ uint16_t PMTWaveformSim::CustomLogNormalPulse(double hit_t0, uint16_t clocktick, double x = (clocktick + fT0Offset) * NS_PER_ADC_SAMPLE - hit_t0; double numerator = log(x/fP1); - numerator = val * val; + numerator = numerator * numerator; double denom = fP2 * fP2; double amplitude = fP0 * exp(-0.5 * numerator/denom) * hit_charge; From 02074e7ea5a73426e7f2e3ff990fdcf8c0c013f9 Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Thu, 27 Feb 2025 11:46:51 -0500 Subject: [PATCH 13/15] Expanding and adjusting BackTracker. The adjustment comes in the use of the MCParticle Idx in all cases rather than the particle ID. This makes access to the underlying particle a bit more straightforward, assuming the vector order does not change. Expansions include the addition of the earliest, mean, and median MCHit times. When Clusters are made with Hits (ie. PMTWaveformSim), there is a matcing step to determine the relevant MCHits, and those are used for these calculations. --- UserTools/BackTracker/BackTracker.cpp | 305 +++++++++++++++++++++----- UserTools/BackTracker/BackTracker.h | 53 ++++- 2 files changed, 300 insertions(+), 58 deletions(-) diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index bf72fc03d..3935b915d 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -27,6 +27,15 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ Log(logmessage, v_error, verbosity); } + bool gotDebugPlots = m_variables.Get("DebugPlots", fDebugPlots); + if (!gotDebugPlots) { + fDebugPlots = false; + logmessage = "BackTracker::Initialize: \"DebugPlots\" not set in the config, defaulting to false"; + Log(logmessage, v_error, verbosity); + } + if (fDebugPlots) + SetupDebug(); + bool gotMCWaveforms = m_variables.Get("MCWaveforms", fMCWaveforms); if (!gotMCWaveforms) { fMCWaveforms = false; @@ -37,7 +46,7 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ // Set up the pointers we're going to save. No need to // delete them at Finalize, the store will handle it - fClusterToBestParticleID = new std::map; + fClusterToBestParticleIdx = new std::map; fClusterToBestParticlePDG = new std::map; fClusterEfficiency = new std::map; fClusterPurity = new std::map; @@ -54,11 +63,14 @@ bool BackTracker::Execute() if (load_status == 2) return true; - fClusterToBestParticleID ->clear(); + fClusterToBestParticleIdx->clear(); fClusterToBestParticlePDG->clear(); fClusterEfficiency ->clear(); fClusterPurity ->clear(); fClusterTotalCharge ->clear(); + fClusterEarliestMCTime ->clear(); + fClusterMeanMCTime ->clear(); + fClusterMedianMCTime ->clear(); fParticleToTankTotalCharge.clear(); @@ -95,51 +107,202 @@ bool BackTracker::Execute() return false; } + // Extract all the MCHits that are in time with this cluster hit 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; + int prtIdx = -5; + int prtPdg = 0; double eff = -5; double pur = -5; double totalCharge = 0; - MatchMCParticle(mcHits, prtId, prtPdg, eff, pur, totalCharge); - - fClusterToBestParticleID ->emplace(apair.first, prtId); + MatchMCParticle(mcHits, prtIdx, prtPdg, eff, pur, totalCharge); + + // Go through the mcHits and determine the earliest, mean, and median hit times + double earliestTime = 99999; + double meanTime = 0; + double medianTime = -5; + for (auto mcHit : mcHits) { + double tempTime = mcHit.GetTime(); + if (tempTime < earliestTime) earliestTime = tempTime; + meanTime += tempTime; + }// end loop over MCHits + if (mcHit.size()) { + meanTime = meanTime / mcHit.size() + + if (mcHit.size() %2 != 0) + medianTime = mcHits.at(mcHit.size()/2); + else + medianTime = (mcHits.at((mcHit.size()-1)/2) + mcHits.at((mcHit.size())/2)) / 2; + }// end mean/median calculation + + fClusterToBestParticleIdx->emplace(apair.first, prtIdx); fClusterToBestParticlePDG->emplace(apair.first, prtPdg); fClusterEfficiency ->emplace(apair.first, eff); fClusterPurity ->emplace(apair.first, pur); fClusterTotalCharge ->emplace(apair.first, totalCharge); - + fClusterEarliestMCTime ->emplace(apair.first, earliestTime); + fClusterMeanMCTime ->emplace(apair.first, meanTime); + fClusterMedianMCTime ->emplace(apair.first, medianTime); + m_data->Stores.at("ANNIEEvent")->Set("MapChannelToPulseTimeToMCHitIdx", fMapChannelToPulseTimeToMCHitIdx ); + + + if (fDebugPlots) { + // Basic cluster values + fDbgClusterTime = apair.first; + fDbgClusterNHits = apair.second.size(); + fDbgClusterNMCHits = mcHits.size(); + fDbgClusterBestParticleIdx = prtIdx; + fDbgClusterBestParticlePDG = prtPdg; + fDbgClusterBestParticleCharge = fParticleToTankTotalCharge.at(prtIdx); + fDbgClusterBestParticleStartEnergy = fMCParticles->at(prtIdx).GetStartEnergy(); + fDbgClusterBestParticleStopEnergy = fMCParticles->at(prtIdx).GetStopEnergy(); + fDbgClusterBestParticleStartTime = fMCParticles->at(prtIdx).GetStartTime(); + fDbgClusterBestParticleStopTime = fMCParticles->at(prtIdx).GetStopTime(); + fDbgClusterEfficiency = eff; + fDbgClusterPurity = pur; + fDbgClusterTotalCharge = totalCharge; + + // Vectors for individual cluster hits + fDbgClusterHitChannel.clear(); + fDbgClusterHitCharge.clear(); + fDbgClusterHitTime.clear(); + for (auto hit : apair.second) { + fDbgClusterHitChannel.push_back(hit.GetTubeId()); + fDbgClusterHitCharge.push_back(hit.GetCharge()); + fDbgClusterHitTime.push_back(hit.GetTime()); + } + + // Vectors for individual matched MCHits and their parent particles + fDbgClusterMCHitChannel.clear(); + fDbgClusterMCHitCharge.clear(); + fDbgClusterMCHitTime.clear(); + fDbgClusterParticleIdx.clear(); + fDbgClusterParticleCharge.clear(); + fDbgClusterParticlePdgCode.clear(); + fDbgClusterParticleStartEnergy.clear(); + fDbgClusterParticleStopEnergy.clear(); + fDbgClusterParticleStartTime.clear(); + fDbgClusterParticleStopTime.clear(); + fDbgClusterMCHitParticleTimeDiff.clear(); + std::vector seenParentIdxs; + for (auto hit : mcHits) { + fDbgClusterMCHitChannel.push_back(hit.GetTubeId()); + fDbgClusterMCHitCharge.push_back(hit.GetCharge()); + fDbgClusterMCHitTime.push_back(hit.GetTime()); + + // Vectors for all matched particles + std::vector parentIdxs = *(hit.GetParents()); + if (!parentIdxs.size()) continue; + for (int parentIdx : parentIdxs) { + if (std::find(seenParentIdxs.begin(), seenParentIdxs.end(), parentIdx) != seenParentIdxs.end()) + continue; // we've recorded info on this particle already + + fDbgClusterParticleIdx.push_back(parentIdx); + fDbgClusterParticleCharge.push_back(fParticleToTankTotalCharge.at(parentIdx)); + fDbgClusterParticlePdgCode.push_back(fMCParticles->at(parentIdx).GetPdgCode()); + fDbgClusterParticleStartEnergy.push_back(fMCParticles->at(parentIdx).GetStartEnergy()); + fDbgClusterParticleStopEnergy.push_back(fMCParticles->at(parentIdx).GetStopEnergy()); + fDbgClusterParticleStartTime.push_back(fMCParticles->at(parentIdx).GetStartTime()); + fDbgClusterParticleStopTime.push_back(fMCParticles->at(parentIdx).GetStopTime()); + + fDbgClusterMCHitParticleTimeDiff.push_back(hit.GetTime() - fMCParticles->at(parentIdx).GetStopTime()); + }// end loop over parent particles + }// end loop over matched MCHits + + + fDbgClusterTree->Fill(); + }// endif debug plots + }// 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 prtIdx = -5; int prtPdg = -5; double eff = -5; double pur = -5; double totalCharge = 0; - MatchMCParticle(apair.second, prtId, prtPdg, eff, pur, totalCharge); + MatchMCParticle(apair.second, prtIdx, prtPdg, eff, pur, totalCharge); - fClusterToBestParticleID ->emplace(apair.first, prtId); + fClusterToBestParticleIdx->emplace(apair.first, prtIdx); fClusterToBestParticlePDG->emplace(apair.first, prtPdg); fClusterEfficiency ->emplace(apair.first, eff); fClusterPurity ->emplace(apair.first, pur); fClusterTotalCharge ->emplace(apair.first, totalCharge); + + if (fDebugPlots) { + // Basic cluster values + fDbgClusterTime = apair.first; + fDbgClusterNMCHits = apair.second.size(); + fDbgClusterBestParticleIdx = prtIdx; + fDbgClusterBestParticlePDG = prtPdg; + fDbgClusterBestParticleCharge = fParticleToTankTotalCharge.at(prtIdx); + fDbgClusterBestParticleStartEnergy = fMCParticles->at(prtIdx).GetStartEnergy(); + fDbgClusterBestParticleStopEnergy = fMCParticles->at(prtIdx).GetStopEnergy(); + fDbgClusterBestParticleStartTime = fMCParticles->at(prtIdx).GetStartTime(); + fDbgClusterBestParticleStopTime = fMCParticles->at(prtIdx).GetStopTime(); + fDbgClusterEfficiency = eff; + fDbgClusterPurity = pur; + fDbgClusterTotalCharge = totalCharge; + + // Vectors for individual matched MCHits and their parent particles + fDbgClusterMCHitChannel.clear(); + fDbgClusterMCHitCharge.clear(); + fDbgClusterMCHitTime.clear(); + fDbgClusterParticleIdx.clear(); + fDbgClusterParticleCharge.clear(); + fDbgClusterParticlePdgCode.clear(); + fDbgClusterParticleStartEnergy.clear(); + fDbgClusterParticleStopEnergy.clear(); + fDbgClusterParticleStartTime.clear(); + fDbgClusterParticleStopTime.clear(); + fDbgClusterMCHitParticleTimeDiff.clear(); + std::vector seenParentIdxs; + for (auto hit : apair.second) { + fDbgClusterMCHitChannel.push_back(hit.GetTubeId()); + fDbgClusterMCHitCharge.push_back(hit.GetCharge()); + fDbgClusterMCHitTime.push_back(hit.GetTime()); + + // Vectors for all matched particles + std::vector parentIdxs = *(hit.GetParents()); + if (!parentIdxs.size()) continue; + for (int parentIdx : parentIdxs) { + if (std::find(seenParentIdxs.begin(), seenParentIdxs.end(), parentIdx) != seenParentIdxs.end()) + continue; // we've recorded info on this particle already + + fDbgClusterParticleIdx.push_back(parentIdx); + fDbgClusterParticleCharge.push_back(fParticleToTankTotalCharge.at(parentIdx)); + fDbgClusterParticlePdgCode.push_back(fMCParticles->at(parentIdx).GetPdgCode()); + fDbgClusterParticleStartEnergy.push_back(fMCParticles->at(parentIdx).GetStartEnergy()); + fDbgClusterParticleStopEnergy.push_back(fMCParticles->at(parentIdx).GetStopEnergy()); + fDbgClusterParticleStartTime.push_back(fMCParticles->at(parentIdx).GetStartTime()); + fDbgClusterParticleStopTime.push_back(fMCParticles->at(parentIdx).GetStopTime()); + + fDbgClusterMCHitParticleTimeDiff.push_back(hit.GetTime() - fMCParticles->at(parentIdx).GetStopTime()); + }// end loop over parent particles + }// end loop over matched MCHits + + + fDbgClusterTree->Fill(); + }// endif debug plots + }// 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("ClusterToBestParticleIdx", fClusterToBestParticleIdx); m_data->Stores.at("ANNIEEvent")->Set("ClusterToBestParticlePDG", fClusterToBestParticlePDG); m_data->Stores.at("ANNIEEvent")->Set("ClusterEfficiency", fClusterEfficiency ); m_data->Stores.at("ANNIEEvent")->Set("ClusterPurity", fClusterPurity ); m_data->Stores.at("ANNIEEvent")->Set("ClusterTotalCharge", fClusterTotalCharge ); + m_data->Stores.at("ANNIEEvent")->Set("ClusterEarliestMCTime", fClusterEarliestMCTime ); + m_data->Stores.at("ANNIEEvent")->Set("ClusterMeanMCTime", fClusterMeanMCTime ); + m_data->Stores.at("ANNIEEvent")->Set("ClusterMedianMCTime", fClusterMedianMCTime ); return true; } @@ -148,6 +311,11 @@ bool BackTracker::Execute() bool BackTracker::Finalise() { + if (fDebugPlots) { + fDebugFile->cd(); + fDbgClusterTree->Write(); + } + return true; } @@ -158,83 +326,70 @@ void BackTracker::SumParticleTankCharge() 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 - // skip any cases we come across + // however, if they do then split the energy equally amongg them std::vector parentIdxs = *(mcHits[mcHitIdx].GetParents()); - if (parentIdxs.size() != 1) continue; - - int particleId = -5; - for (auto bpair : *fMCParticleIndexMap) { - if (bpair.second == parentIdxs[0]) particleId = bpair.first; - } - if (particleId == -5) continue; - - double depositedCharge = mcHits[mcHitIdx].GetCharge(); - if (!fParticleToTankTotalCharge.count(particleId)) - fParticleToTankTotalCharge.emplace(particleId, depositedCharge); + if (!parentIdxs.size()) continue; // no parents recorded? + + double depositedCharge = mcHits[mcHitIdx].GetCharge() / parentIdxs.size(); + for (int parentIdx : parentIdxs) { + if (!fParticleToTankTotalCharge.count(parentIdx)) + fParticleToTankTotalCharge.emplace(parentIdx, depositedCharge); else - fParticleToTankTotalCharge.at(particleId) += depositedCharge; - } - } + fParticleToTankTotalCharge.at(parentIdx) += depositedCharge; + }// end loop over parent indexes + }// end loop over MCHits + }// end loop over PMTs } //------------------------------------------------------------------------------ -void BackTracker::MatchMCParticle(std::vector const &mchits, int &prtId, int &prtPdg, double &eff, double &pur, double &totalCharge) +void BackTracker::MatchMCParticle(std::vector const &mchits, int &prtIdx, int &prtPdg, double &eff, double &pur, double &totalCharge) { // Loop over the hits and get all of their parents and the energy that each one contributed - // be sure to bunch up all neutronic contributions std::map mapParticleToTotalClusterCharge; totalCharge = 0; for (auto mchit : mchits) { std::vector parentIdxs = *(mchit.GetParents()); - if (parentIdxs.size() != 1) { - logmessage = "BackTracker::MatchMCParticle: this MCHit has "; - logmessage += std::to_string(parentIdxs.size()) + " parents!"; - Log(logmessage, v_debug, verbosity); - continue; - } - - int particleId = -5; - for (auto apair : *fMCParticleIndexMap) { - if (apair.second == parentIdxs[0]) particleId = apair.first; - } - if (particleId == -5) continue; - - double depositedCharge = mchit.GetCharge(); - totalCharge += depositedCharge; + if (!parentIdxs.size()) continue; + + double hitCharge = mchit.GetCharge(); + totalCharge += hitCharge; + hitCharge = hitCharge / parentIdxs.size(); - if (mapParticleToTotalClusterCharge.count(particleId) == 0) - mapParticleToTotalClusterCharge.emplace(particleId, depositedCharge); - else - mapParticleToTotalClusterCharge[particleId] += depositedCharge; - } + for (int parentIdx : parentIdxs) { + if (mapParticleToTotalClusterCharge.count(parentIdx) == 0) + mapParticleToTotalClusterCharge.emplace(parentIdx, hitCharge); + else + mapParticleToTotalClusterCharge[parentIdx] += hitCharge; + }// end loop over parentIdxs + }// end loop over MCHits // Loop over the particleIds to find the primary contributer to the cluster double maxCharge = 0; for (auto apair : mapParticleToTotalClusterCharge) { if (apair.second > maxCharge) { maxCharge = apair.second; - prtId = apair.first; + prtIdx = apair.first; } - } + }// end loop over the particle charge map // Check that we have some charge, if not then something is wrong so pass back all -5 if (totalCharge > 0) { - eff = maxCharge/fParticleToTankTotalCharge.at(prtId); + eff = maxCharge/fParticleToTankTotalCharge.at(prtIdx); pur = maxCharge/totalCharge; - prtPdg = (fMCParticles->at(fMCParticleIndexMap->at(prtId))).GetPdgCode(); + prtPdg = (fMCParticles->at(prtIdx)).GetPdgCode(); } else { - prtId = -5; + prtIdx = -5; eff = -5; pur = -5; totalCharge = -5; } logmessage = "BackTracker::MatchMCParticle: best particleId is : "; - logmessage += std::to_string(prtId) + " which has PDG: " + std::to_string(prtPdg); + logmessage += std::to_string(prtIdx) + " which has PDG: " + std::to_string(prtPdg); Log(logmessage, v_message, verbosity); - } //------------------------------------------------------------------------------ @@ -274,6 +429,7 @@ bool BackTracker::MapPulsesToParentIdxs() 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()) @@ -296,6 +452,43 @@ bool BackTracker::MapPulsesToParentIdxs() return true; } +//------------------------------------------------------------------------------ +void BackTracker::SetupDebug() +{ + fDebugFile = new TFile("BackTracker_Debug.root", "RECREATE"); + fDbgClusterTree = new TTree("clusterTree", "clusterTree"); + + fDbgClusterTree->Branch("Time", &fDbgClusterTime); + fDbgClusterTree->Branch("NHits", &fDbgClusterNHits); + fDbgClusterTree->Branch("NMCHits", &fDbgClusterNMCHits); + fDbgClusterTree->Branch("BestParticleIdx", &fDbgClusterBestParticleIdx); + fDbgClusterTree->Branch("BestParticlePDG", &fDbgClusterBestParticlePDG); + fDbgClusterTree->Branch("BestParticleCharge", &fDbgClusterBestParticleCharge); + fDbgClusterTree->Branch("BestParticleStartEnergy", &fDbgClusterBestParticleStartEnergy); + fDbgClusterTree->Branch("BestParticleStopEnergy", &fDbgClusterBestParticleStopEnergy); + fDbgClusterTree->Branch("BestParticleStartTime", &fDbgClusterBestParticleStartTime); + fDbgClusterTree->Branch("BestParticleStopTime", &fDbgClusterBestParticleStopTime); + fDbgClusterTree->Branch("Efficiency", &fDbgClusterEfficiency); + fDbgClusterTree->Branch("Purity", &fDbgClusterPurity); + fDbgClusterTree->Branch("TotalCharge", &fDbgClusterTotalCharge); + fDbgClusterTree->Branch("HitChannel", &fDbgClusterHitChannel); + fDbgClusterTree->Branch("HitCharge", &fDbgClusterHitCharge); + fDbgClusterTree->Branch("HitTime", &fDbgClusterHitTime); + fDbgClusterTree->Branch("MCHitChannel", &fDbgClusterMCHitChannel); + fDbgClusterTree->Branch("MCHitCharge", &fDbgClusterMCHitCharge); + fDbgClusterTree->Branch("MCHitTime", &fDbgClusterMCHitTime); + fDbgClusterTree->Branch("ParticleIdx", &fDbgClusterParticleIdx); + fDbgClusterTree->Branch("ParticleCharge", &fDbgClusterParticleCharge); + fDbgClusterTree->Branch("ParticlePDG", &fDbgClusterParticlePdgCode); + fDbgClusterTree->Branch("ParticleStartEnergy", &fDbgClusterParticleStartEnergy); + fDbgClusterTree->Branch("ParticleStopEnergy", &fDbgClusterParticleStopEnergy); + fDbgClusterTree->Branch("ParticleStartTime", &fDbgClusterParticleStartTime); + fDbgClusterTree->Branch("ParticleStopTime", &fDbgClusterParticleStopTime); + fDbgClusterTree->Branch("MCHitParticleTimeDiff", &fDbgClusterMCHitParticleTimeDiff); + +} + + //------------------------------------------------------------------------------ int BackTracker::LoadFromStores() { diff --git a/UserTools/BackTracker/BackTracker.h b/UserTools/BackTracker/BackTracker.h index 31566dc47..56989d110 100644 --- a/UserTools/BackTracker/BackTracker.h +++ b/UserTools/BackTracker/BackTracker.h @@ -9,6 +9,10 @@ #include "Particle.h" #include "ADCPulse.h" +#include "TFile.h" +#include "TTree.h" + + /** * \class BackTracker @@ -31,9 +35,12 @@ class BackTracker: public Tool { int LoadFromStores(); ///< Does all the loading so I can move it away from the Execute function void SumParticleTankCharge(); - void MatchMCParticle(std::vector const &mchits, int &prtId, int &prtPdg, double &eff, double &pur, double &totalCharge); ///< The meat and potatoes + void MatchMCParticle(std::vector const &mchits, int &prtIdx, int &prtPdg, double &eff, double &pur, double &totalCharge); ///< The meat and potatoes bool MapPulsesToParentIdxs(); + + void SetupDebug(); + private: @@ -62,13 +69,55 @@ class BackTracker: public Tool { // the the purity based on the best matched particle // the total deposited charge in the cluster // the ammount of cluster charge due to neutrons - std::map *fClusterToBestParticleID = nullptr; + std::map *fClusterToBestParticleIdx = nullptr; std::map *fClusterToBestParticlePDG = nullptr; std::map *fClusterEfficiency = nullptr; std::map *fClusterPurity = nullptr; std::map *fClusterTotalCharge = nullptr; + std::map *fClusterEarliestMCTime = nullptr; + std::map *fClusterMeanMCTime = nullptr; + std::map *fClusterMedianMCTime = nullptr; + // Config variables + bool fDebugPlots; bool fMCWaveforms; + + // For debug plots + TFile *fDebugFile; + + TTree *fDbgClusterTree; + double fDbgClusterTime; + int fDbgClusterNHits; + int fDbgClusterNMCHits; + int fDbgClusterBestParticleIdx; + int fDbgClusterBestParticlePDG; + double fDbgClusterBestParticleCharge; + double fDbgClusterBestParticleStartEnergy; + double fDbgClusterBestParticleStopEnergy; + double fDbgClusterBestParticleStartTime; + double fDbgClusterBestParticleStopTime; + double fDbgClusterEfficiency; + double fDbgClusterPurity; + double fDbgClusterTotalCharge; + std::vector fDbgClusterHitChannel; + std::vector fDbgClusterHitCharge; + std::vector fDbgClusterHitTime; + std::vector fDbgClusterMCHitChannel; + std::vector fDbgClusterMCHitCharge; + std::vector fDbgClusterMCHitTime; + std::vector fDbgClusterParticleIdx; + std::vector fDbgClusterParticleCharge; + std::vector fDbgClusterParticlePdgCode; + std::vector fDbgClusterParticleStartEnergy; + std::vector fDbgClusterParticleStopEnergy; + std::vector fDbgClusterParticleStartTime; + std::vector fDbgClusterParticleStopTime; + std::vector fDbgClusterMCHitParticleTimeDiff; + + + + + /// \brief verbosity levels: if 'verbosity' < this level, the message type will be logged. int verbosity; From daed73210ac23c4a50045503050930fe0152f795 Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Thu, 27 Feb 2025 14:00:32 -0500 Subject: [PATCH 14/15] Looks like I pushed a buggy version of this. Fixing. --- UserTools/BackTracker/BackTracker.cpp | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index 3935b915d..863da7c0a 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -130,13 +130,16 @@ bool BackTracker::Execute() if (tempTime < earliestTime) earliestTime = tempTime; meanTime += tempTime; }// end loop over MCHits - if (mcHit.size()) { - meanTime = meanTime / mcHit.size() + if (mcHits.size()) { + meanTime = meanTime / mcHits.size(); - if (mcHit.size() %2 != 0) - medianTime = mcHits.at(mcHit.size()/2); - else - medianTime = (mcHits.at((mcHit.size()-1)/2) + mcHits.at((mcHit.size())/2)) / 2; + if (mcHits.size() %2 != 0) + medianTime = mcHits.at(mcHits.size()/2).GetTime(); + else { + double time1 = mcHits.at((mcHits.size()-1)/2).GetTime(); + double time2 = mcHits.at((mcHits.size())/2).GetTime(); + medianTime = (time1 + time2) / 2; + } }// end mean/median calculation fClusterToBestParticleIdx->emplace(apair.first, prtIdx); From aeba7bb3867c7b7ae4a188aaa342108494f6dbdb Mon Sep 17 00:00:00 2001 From: Andrew Sutton Date: Thu, 27 Feb 2025 15:37:44 -0500 Subject: [PATCH 15/15] Forgot to initialize the new pointers... --- UserTools/BackTracker/BackTracker.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index 863da7c0a..57830c390 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -51,6 +51,10 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ fClusterEfficiency = new std::map; fClusterPurity = new std::map; fClusterTotalCharge = new std::map; + fClusterEarliestMCTime = new std::map; + fClusterMeanMCTime = new std::map; + fClusterMedianMCTime = new std::map; + return true; } @@ -68,9 +72,9 @@ bool BackTracker::Execute() fClusterEfficiency ->clear(); fClusterPurity ->clear(); fClusterTotalCharge ->clear(); - fClusterEarliestMCTime ->clear(); - fClusterMeanMCTime ->clear(); - fClusterMedianMCTime ->clear(); + fClusterEarliestMCTime ->clear(); + fClusterMeanMCTime ->clear(); + fClusterMedianMCTime ->clear(); fParticleToTankTotalCharge.clear();