From 9c965801ea13d1ddada96bb61233cc1770d0546f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Mon, 6 Jul 2026 08:59:36 -0500 Subject: [PATCH 1/5] Modify light propagation correction to include particle propagation time --- .../LightPropagationCorrectionAna_module.cc | 10 +- .../LightPropagationCorrection_module.cc | 115 ++++++++++++++++-- .../LightPropagationCorrection_module.hh | 13 +- ...sbnd_lightpropagationcorrection_config.fcl | 4 + 4 files changed, 130 insertions(+), 12 deletions(-) diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc b/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc index 9b09202a2..2f29e3662 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc @@ -56,6 +56,8 @@ class LightPropagationCorrectionAna : public art::EDAnalyzer { double fNuToFLight; double fNuToFCharge; double fOpFlashT0Corrected; + double fParticlePropagationTime; + double fPhotonPropagationTime; unsigned int _eventID; unsigned int _runID; @@ -95,7 +97,6 @@ void LightPropagationCorrectionAna::analyze(art::Event const& e) return; } - std::vector> correctedopflash_v; art::fill_ptr_vector(correctedopflash_v, correctedopflash_h); @@ -121,11 +122,15 @@ void LightPropagationCorrectionAna::analyze(art::Event const& e) std::cout << "Corrected flash time light only " << correctedopflash->NuToFLight << std::endl; std::cout << "Corrected flash time tpc z corr " << correctedopflash->NuToFCharge << std::endl; std::cout << "Corrected flash time prop corr tpc z corr " << correctedopflash->OpFlashT0Corrected << std::endl; + std::cout << "Particle propagation time " << correctedopflash->ParticlePropagationTime << std::endl; + std::cout << "Photon propagation time " << correctedopflash->PhotonPropagationTime << std::endl; fOpFlashT0 = correctedopflash->OpFlashT0; fNuToFLight = correctedopflash->NuToFLight; fNuToFCharge = correctedopflash->NuToFCharge; fOpFlashT0Corrected = correctedopflash->OpFlashT0Corrected; + fParticlePropagationTime = correctedopflash->ParticlePropagationTime; + fPhotonPropagationTime = correctedopflash->PhotonPropagationTime; fTree->Fill(); } } @@ -144,7 +149,8 @@ void LightPropagationCorrectionAna::beginJob() fTree->Branch("fNuToFLight", &fNuToFLight, "NuToFLight/d"); fTree->Branch("fNuToFCharge", &fNuToFCharge, "NuToFCharge/d"); fTree->Branch("fOpFlashT0Corrected", &fOpFlashT0Corrected, "OpFlashT0Corrected/d"); - + fTree->Branch("fParticlePropagationTime", &fParticlePropagationTime, "ParticlePropagationTime/d"); + fTree->Branch("fPhotonPropagationTime", &fPhotonPropagationTime, "PhotonPropagationTime/d"); } DEFINE_ART_MODULE(LightPropagationCorrectionAna) \ No newline at end of file diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc index 6b3e999bc..b32ebd163 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc @@ -22,6 +22,9 @@ sbnd::LightPropagationCorrection::LightPropagationCorrection(fhicl::ParameterSet fVGroupVUV( p.get("VGroupVUV") ), fNuScoreThreshold( p.get("NuScoreThreshold") ), fFMScoreThreshold( p.get("FMScoreThreshold") ), + fMinHitPE ( p.get("MinHitPE") ), + fPreWindow ( p.get("PreWindow") ), + fPostWindow ( p.get("PostWindow") ), fDebug( p.get("Debug", false) ) // // More initializers here. @@ -29,6 +32,8 @@ sbnd::LightPropagationCorrection::LightPropagationCorrection(fhicl::ParameterSet fNOpChannels = fWireReadout.NOpChannels(); // Initialize the TimeCorrectionVector PerChannel fTimeCorrectionPerChannel.resize(fNOpChannels, 0.0); // Initialize with zero or any default value + fParticlePropagationTimePerChannel.resize(fNOpChannels, 0.0); // Initialize with zero or any default value + fPhotonPropagationTimePerChannel.resize(fNOpChannels, 0.0); // Initialize with zero or any default value for(unsigned int opch=0; opchTime(); double newFlashTime = 0.0; + double particlePropTime = 0.0; + double photonPropTime = 0.0; for(const auto& lflash : flash_v) { // Get Flash Barycenter double Ycenter, Zcenter, Ywidth, Zwidth; @@ -267,16 +274,21 @@ void sbnd::LightPropagationCorrection::produce(art::Event & e) double flasht0 = lflash.time; // Refine t0 calculation flasht0 = _flasht0calculator->GetFlashT0(lflash.time, GetAssociatedLiteHits(lflash, ophits)); + this->GetSelectedChannelsFlash(lflash.time, GetAssociatedLiteHits(lflash, ophits)); recob::OpFlash flash(flasht0, lflash.time_err, flasht0, ( flasht0) / 1600., lflash.channel_pe, 0, 0, 1, // this are just default values 100., -1., Ycenter, Ywidth, Zcenter, Zwidth); newFlashTime = flasht0; + particlePropTime = GetAverageParticlePropagationTime()/1000; + photonPropTime = GetAveragePhotonPropagationTime()/1000; sbn::CorrectedOpFlashTiming correctedOpFlashTiming; correctedOpFlashTiming.OpFlashT0 = originalFlashTime; correctedOpFlashTiming.NuToFLight = (Zcenter/fSpeedOfLight)/1000; correctedOpFlashTiming.NuToFCharge = (fRecoVz/fSpeedOfLight)/1000; correctedOpFlashTiming.OpFlashT0Corrected = newFlashTime; + correctedOpFlashTiming.ParticlePropagationTime = particlePropTime; + correctedOpFlashTiming.PhotonPropagationTime = photonPropTime; correctedOpFlashTimes->emplace_back(std::move(correctedOpFlashTiming)); } @@ -396,6 +408,8 @@ void sbnd::LightPropagationCorrection::ResetSliceInfo() fSpacePointZ.clear(); fSpacePointIntegral.clear(); fTimeCorrectionPerChannel.assign(fNOpChannels, 0.0); // Reset the time correction vector for each channel + fParticlePropagationTimePerChannel.assign(fNOpChannels, 0.0); // Reset the particle propagation time vector for each channel + fPhotonPropagationTimePerChannel.assign(fNOpChannels, 0.0); // Reset the photon propagation time vector for each channel fChargeBarycenterX.assign(2, 0.0); fChargeBarycenterY.assign(2, 0.0); fChargeBarycenterZ.assign(2, 0.0); @@ -408,13 +422,17 @@ void sbnd::LightPropagationCorrection::ResetSliceInfo() void sbnd::LightPropagationCorrection::GetPropagationTimeCorrectionPerChannel() { // Implementation - for(size_t opdet = 0; opdet < fOpDetID.size(); ++opdet) { +for(size_t opdet = 0; opdet < fOpDetID.size(); ++opdet) { double _opDetX = fOpDetX[opdet]; double _opDetY = fOpDetY[opdet]; double _opDetZ = fOpDetZ[opdet]; float minPropTime = 999999999.; + float minPartPropTime = 999999999.; + float minLightPropTime = 999999999.; for(size_t sp=0; sp 0; + if(!isInSameTPC) continue; // Skip points not in the same TPC double dx = fSpacePointX[sp] - _opDetX; double dy = fSpacePointY[sp] - _opDetY; double dz = fSpacePointZ[sp] - _opDetZ; @@ -425,18 +443,21 @@ void sbnd::LightPropagationCorrection::GetPropagationTimeCorrectionPerChannel() float lightPropTimeVIS = spToCathode/fVGroupVUV + cathodeToOpDet/fVGroupVIS; // Speed float lightPropTimeVUV = distanceToOpDet / fVGroupVUV; // Speed of light in mm/ns for VUV float lightPropTime = 0; - const std::string pdType = fPDSMap.pdType(opdet); - if(pdType=="pmt_coated" || pdType=="xarapuca_vuv") + if(fPDSMap.pdType(opdet)=="pmt_coated" || fPDSMap.pdType(opdet)=="xarapuca_vuv") lightPropTime = std::min(lightPropTimeVIS, lightPropTimeVUV); - else if(pdType=="pmt_uncoated" || pdType=="xarapuca_vis") + else if(fPDSMap.pdType(opdet)=="pmt_uncoated" || fPDSMap.pdType(opdet)=="xarapuca_vis") lightPropTime = lightPropTimeVIS; - else - throw std::runtime_error("LightPropagationCorrection: unexpected pdType '" + pdType + "' for opdet " + std::to_string(opdet)); float partPropTime = std::sqrt((fSpacePointX[sp]-fRecoVx)*(fSpacePointX[sp]-fRecoVx) + (fSpacePointY[sp]-fRecoVy)*(fSpacePointY[sp]-fRecoVy) + (fSpacePointZ[sp]-fRecoVz)*(fSpacePointZ[sp]-fRecoVz))/fSpeedOfLight; float PropTime = lightPropTime + partPropTime; - if(PropTime < minPropTime) minPropTime = PropTime; + if(PropTime < minPropTime) { + minPropTime = PropTime; + minPartPropTime = partPropTime; + minLightPropTime = lightPropTime; + } } fTimeCorrectionPerChannel[opdet] = -minPropTime; + fParticlePropagationTimePerChannel[opdet] = minPartPropTime; + fPhotonPropagationTimePerChannel[opdet] = minLightPropTime; } } @@ -515,7 +536,7 @@ void sbnd::LightPropagationCorrection::FillCorrectionTree(double & newFlashTime, ::lightana::LiteOpHitArray_t sbnd::LightPropagationCorrection::GetAssociatedLiteHits(::lightana::LiteOpFlash_t lite_flash, ::lightana::LiteOpHitArray_t lite_hits_v) - { +{ ::lightana::LiteOpHitArray_t flash_hits_v; for(auto const& hitidx : lite_flash.asshit_idx) { @@ -523,4 +544,80 @@ ::lightana::LiteOpHitArray_t sbnd::LightPropagationCorrection::GetAssociatedLite } return flash_hits_v; - } +} + +double sbnd::LightPropagationCorrection::GetAverageParticlePropagationTime() +{ + double sum = 0.0; + int n = 0; + for (size_t i=0; i > selected_hits; + double pe_sum = 0.0; + + // Limpiar vector de salida por si ya tenía contenido + fSelectedChannelList.clear(); + // Fill vector with selected hits in the specified window + for (auto const& hit : ophit_list) { + + if (hit.peak_time < flash_time + fPostWindow && + hit.peak_time > flash_time - fPreWindow && + hit.pe > fMinHitPE) { + + // (PE, peak_time, channel) + selected_hits.emplace_back(hit.pe, hit.peak_time, hit.channel); + pe_sum += hit.pe; + } + } + + // Sort by PE in descending order + std::sort( + selected_hits.begin(), + selected_hits.end(), + [](auto const& a, auto const& b) { + return std::get<0>(a) > std::get<0>(b); + } + ); + + double pe_count = 0.0; + + // Loop over selected ophits + for (size_t ix = 0; ix < selected_hits.size(); ix++) { + + double pe = std::get<0>(selected_hits[ix]); + size_t channel = std::get<2>(selected_hits[ix]); + + // Guardar channel ID + fSelectedChannelList.push_back(channel); + + pe_count += pe; + if (pe_count / pe_sum > .5) + break; + } +} \ No newline at end of file diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh index d0ae118fd..dc33c9ef3 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh @@ -142,7 +142,9 @@ private: void FillLiteOpHit(std::vector const& , std::vector<::lightana::LiteOpHit_t>& ); void FillCorrectionTree(double & , recob::OpFlash const& , std::vector const& , std::vector const& ); ::lightana::LiteOpHitArray_t GetAssociatedLiteHits(::lightana::LiteOpFlash_t , ::lightana::LiteOpHitArray_t ); - + void GetSelectedChannelsFlash(double , ::lightana::LiteOpHitArray_t ); + double GetAverageParticlePropagationTime(); + double GetAveragePhotonPropagationTime(); geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); opdet::sbndPDMapAlg fPDSMap; @@ -176,6 +178,9 @@ private: bool fSaveCorrectionTree; std::vector fTimeCorrectionPerChannel; + std::vector fParticlePropagationTimePerChannel; + std::vector fPhotonPropagationTimePerChannel; + std::vector fSelectedChannelList; double fRecoVx = 0.0; double fRecoVy = 0.0; double fRecoVz = 0.0; @@ -206,6 +211,12 @@ private: double fNuScoreThreshold; double fFMScoreThreshold; + // Parameters for the selected channel computation + int fMinHitPE; + double fPreWindow; + double fPostWindow; + + bool fDebug; art::ServiceHandle tfs; diff --git a/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl index 3277a7bb8..2205072a7 100644 --- a/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl +++ b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl @@ -35,6 +35,10 @@ LightPropagationCorrection: NuScoreThreshold: 0.1 FMScoreThreshold: 0.001 + MinHitPE: @local::FlashT0SelectedChannels.MinHitPE + PreWindow: @local::FlashT0SelectedChannels.PreWindow + PostWindow: @local::FlashT0SelectedChannels.PostWindow + FlashMatchingTool: "BarycenterFM" } From c806e705919b68537fdbee629baa91b99642f3fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Mon, 6 Jul 2026 09:05:07 -0500 Subject: [PATCH 2/5] Add pdfraction as parameter --- .../LightPropagationCorrection_module.cc | 3 ++- .../LightPropagationCorrection_module.hh | 1 + .../job/sbnd_lightpropagationcorrection_config.fcl | 1 + 3 files changed, 4 insertions(+), 1 deletion(-) diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc index b32ebd163..eafe0ee00 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc @@ -25,6 +25,7 @@ sbnd::LightPropagationCorrection::LightPropagationCorrection(fhicl::ParameterSet fMinHitPE ( p.get("MinHitPE") ), fPreWindow ( p.get("PreWindow") ), fPostWindow ( p.get("PostWindow") ), + fPDFraction ( p.get("PDFraction") ), fDebug( p.get("Debug", false) ) // // More initializers here. @@ -617,7 +618,7 @@ void sbnd::LightPropagationCorrection::GetSelectedChannelsFlash( fSelectedChannelList.push_back(channel); pe_count += pe; - if (pe_count / pe_sum > .5) + if (pe_count / pe_sum > fPDFraction) break; } } \ No newline at end of file diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh index dc33c9ef3..1bbe62712 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh @@ -215,6 +215,7 @@ private: int fMinHitPE; double fPreWindow; double fPostWindow; + double fPDFraction; bool fDebug; diff --git a/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl index 2205072a7..ecfd6e2e8 100644 --- a/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl +++ b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl @@ -38,6 +38,7 @@ LightPropagationCorrection: MinHitPE: @local::FlashT0SelectedChannels.MinHitPE PreWindow: @local::FlashT0SelectedChannels.PreWindow PostWindow: @local::FlashT0SelectedChannels.PostWindow + PDFraction: @local::FlashT0SelectedChannels.PDFraction FlashMatchingTool: "BarycenterFM" } From 3e702ac1ec8e630765c4d5109ada921e6c82c0f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Tue, 7 Jul 2026 04:11:57 -0500 Subject: [PATCH 3/5] Clean code adding correctopflash function --- .../LightPropagationCorrection_module.cc | 141 +++++++++--------- .../LightPropagationCorrection_module.hh | 5 +- 2 files changed, 71 insertions(+), 75 deletions(-) diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc index eafe0ee00..e064f715f 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc @@ -119,8 +119,8 @@ void sbnd::LightPropagationCorrection::produce(art::Event & e) //PFP to space points art::FindManyP pfp_sp_assns(pfpHandle, e, fSpacePointLabel); //OpFlash to OpHit - art::FindManyP flashToOpHitAssns_tpc0(opflashListHandle_tpc0, e, fOpFlashLabel_tpc0); - art::FindManyP flashToOpHitAssns_tpc1(opflashListHandle_tpc1, e, fOpFlashLabel_tpc1); + flashToOpHitAssns_tpc0 = std::make_unique>( opflashListHandle_tpc0, e, fOpFlashLabel_tpc0); + flashToOpHitAssns_tpc1 = std::make_unique>(opflashListHandle_tpc1, e, fOpFlashLabel_tpc1); // PFP Metadata art::FindManyP pfp_to_metadata(pfpHandle, e, fReco2Label); @@ -237,70 +237,14 @@ void sbnd::LightPropagationCorrection::produce(art::Event & e) } else throw art::Exception(art::errors::LogicError) << " Flash matching tool " << fFlashMatchingTool << " not supported ." << std::endl; - - // Get the ophits associated to the flash - std::vector> ophitlist; - if(flashFM[0]->XCenter()<0) - { - ophitlist = flashToOpHitAssns_tpc0.at(flashFM[0].key()); - _mgr = _mgr_tpc0; // Use the TPC 0 flash finder manager - } - else - { - ophitlist = flashToOpHitAssns_tpc1.at(flashFM[0].key()); - _mgr = _mgr_tpc1; // Use the TPC 1 flash finder manager - } - - std::vector newOpHitList; - std::vector oldOpHitList; - for(const auto& ophit : ophitlist) { - oldOpHitList.push_back(*ophit); - } - // Get the list of the corrected OpHits - this->CorrectOpHitTime(ophitlist, newOpHitList); - // Create the list of ophit lite to be used in the flash finder - ::lightana::LiteOpHitArray_t ophits; - this->FillLiteOpHit(newOpHitList, ophits); - // Create the flash manager - auto const flash_v = _mgr.RecoFlash(ophits); - double originalFlashTime = flashFM[0]->Time(); - double newFlashTime = 0.0; - double particlePropTime = 0.0; - double photonPropTime = 0.0; - for(const auto& lflash : flash_v) { - // Get Flash Barycenter - double Ycenter, Zcenter, Ywidth, Zwidth; - _flashgeo->GetFlashLocation(lflash.channel_pe, Ycenter, Zcenter, Ywidth, Zwidth); - // Get flasht0 - double flasht0 = lflash.time; - // Refine t0 calculation - flasht0 = _flasht0calculator->GetFlashT0(lflash.time, GetAssociatedLiteHits(lflash, ophits)); - this->GetSelectedChannelsFlash(lflash.time, GetAssociatedLiteHits(lflash, ophits)); - recob::OpFlash flash(flasht0, lflash.time_err, flasht0, - ( flasht0) / 1600., lflash.channel_pe, - 0, 0, 1, // this are just default values - 100., -1., Ycenter, Ywidth, Zcenter, Zwidth); - newFlashTime = flasht0; - particlePropTime = GetAverageParticlePropagationTime()/1000; - photonPropTime = GetAveragePhotonPropagationTime()/1000; - sbn::CorrectedOpFlashTiming correctedOpFlashTiming; - correctedOpFlashTiming.OpFlashT0 = originalFlashTime; - correctedOpFlashTiming.NuToFLight = (Zcenter/fSpeedOfLight)/1000; - correctedOpFlashTiming.NuToFCharge = (fRecoVz/fSpeedOfLight)/1000; - correctedOpFlashTiming.OpFlashT0Corrected = newFlashTime; - correctedOpFlashTiming.ParticlePropagationTime = particlePropTime; - correctedOpFlashTiming.PhotonPropagationTime = photonPropTime; - correctedOpFlashTimes->emplace_back(std::move(correctedOpFlashTiming)); - } - + sbn::CorrectedOpFlashTiming correctedOpFlashTiming; + CorrectOpFlash(flashFM[0], correctedOpFlashTiming); + correctedOpFlashTimes->emplace_back(std::move(correctedOpFlashTiming)); art::Ptr newCorrectedOpFlashTimingPtr = make_correctedopflashtime_ptr(correctedOpFlashTimes->size()-1); newCorrectedOpFlashTimingSliceAssn->addSingle(slice, newCorrectedOpFlashTimingPtr); newCorrectedOpFlashTimingOpFlashAssn->addSingle(flashFM[0], newCorrectedOpFlashTimingPtr); - if(fSaveCorrectionTree){ - this->FillCorrectionTree(newFlashTime, *flashFM[0], oldOpHitList, newOpHitList); - } } - if(fSaveCorrectionTree) fTree->Fill(); + ResetEventVars(); e.put(std::move(correctedOpFlashTimes)); e.put(std::move(newCorrectedOpFlashTimingSliceAssn)); @@ -578,11 +522,8 @@ void sbnd::LightPropagationCorrection::GetSelectedChannelsFlash( double flash_time, ::lightana::LiteOpHitArray_t ophit_list) { - // Ahora guardamos también el channel ID std::vector< std::tuple > selected_hits; double pe_sum = 0.0; - - // Limpiar vector de salida por si ya tenía contenido fSelectedChannelList.clear(); // Fill vector with selected hits in the specified window for (auto const& hit : ophit_list) { @@ -598,14 +539,7 @@ void sbnd::LightPropagationCorrection::GetSelectedChannelsFlash( } // Sort by PE in descending order - std::sort( - selected_hits.begin(), - selected_hits.end(), - [](auto const& a, auto const& b) { - return std::get<0>(a) > std::get<0>(b); - } - ); - + std::sort(selected_hits.begin(), selected_hits.end(), [](auto const& a, auto const& b) {return std::get<0>(a) > std::get<0>(b);}); double pe_count = 0.0; // Loop over selected ophits @@ -614,11 +548,70 @@ void sbnd::LightPropagationCorrection::GetSelectedChannelsFlash( double pe = std::get<0>(selected_hits[ix]); size_t channel = std::get<2>(selected_hits[ix]); - // Guardar channel ID + // Save channel ID fSelectedChannelList.push_back(channel); pe_count += pe; if (pe_count / pe_sum > fPDFraction) break; } +} + +void sbnd::LightPropagationCorrection::CorrectOpFlash(art::Ptr const& flash, sbn::CorrectedOpFlashTiming &correctedOpFlashTiming) +{ + // Get the ophits associated to the flash + std::vector> ophitlist; + if(flash->XCenter()<0) + { + ophitlist = flashToOpHitAssns_tpc0->at(flash.key()); + _mgr = _mgr_tpc0; // Use the TPC 0 flash finder manager + } + else + { + ophitlist = flashToOpHitAssns_tpc1->at(flash.key()); + _mgr = _mgr_tpc1; // Use the TPC 1 flash finder manager + } + std::vector newOpHitList; + std::vector oldOpHitList; + for(const auto& ophit : ophitlist) { + oldOpHitList.push_back(*ophit); + } + // Get the list of the corrected OpHits + this->CorrectOpHitTime(ophitlist, newOpHitList); + // Create the list of ophit lite to be used in the flash finder + ::lightana::LiteOpHitArray_t ophits; + this->FillLiteOpHit(newOpHitList, ophits); + // Create the flash manager + auto const flash_v = _mgr.RecoFlash(ophits); + double originalFlashTime = flash->Time(); + double newFlashTime = 0.0; + double particlePropTime = 0.0; + double photonPropTime = 0.0; + for(const auto& lflash : flash_v) { + // Get Flash Barycenter + double Ycenter, Zcenter, Ywidth, Zwidth; + _flashgeo->GetFlashLocation(lflash.channel_pe, Ycenter, Zcenter, Ywidth, Zwidth); + // Get flasht0 + double flasht0 = lflash.time; + // Refine t0 calculation + flasht0 = _flasht0calculator->GetFlashT0(lflash.time, GetAssociatedLiteHits(lflash, ophits)); + this->GetSelectedChannelsFlash(lflash.time, GetAssociatedLiteHits(lflash, ophits)); + recob::OpFlash flash(flasht0, lflash.time_err, flasht0, + ( flasht0) / 1600., lflash.channel_pe, + 0, 0, 1, // this are just default values + 100., -1., Ycenter, Ywidth, Zcenter, Zwidth); + newFlashTime = flasht0; + particlePropTime = GetAverageParticlePropagationTime()/1000; + photonPropTime = GetAveragePhotonPropagationTime()/1000; + correctedOpFlashTiming.OpFlashT0 = originalFlashTime; + correctedOpFlashTiming.NuToFLight = (Zcenter/fSpeedOfLight)/1000; + correctedOpFlashTiming.NuToFCharge = (fRecoVz/fSpeedOfLight)/1000; + correctedOpFlashTiming.OpFlashT0Corrected = newFlashTime; + correctedOpFlashTiming.ParticlePropagationTime = particlePropTime; + correctedOpFlashTiming.PhotonPropagationTime = photonPropTime; + } + + if(fSaveCorrectionTree){ + this->FillCorrectionTree(newFlashTime, *flash, oldOpHitList, newOpHitList); + } } \ No newline at end of file diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh index 1bbe62712..55b5944ff 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh @@ -145,9 +145,13 @@ private: void GetSelectedChannelsFlash(double , ::lightana::LiteOpHitArray_t ); double GetAverageParticlePropagationTime(); double GetAveragePhotonPropagationTime(); + void CorrectOpFlash(art::Ptr const& , sbn::CorrectedOpFlashTiming &); geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); opdet::sbndPDMapAlg fPDSMap; + // Flash ophit associations + std::unique_ptr> flashToOpHitAssns_tpc0; + std::unique_ptr> flashToOpHitAssns_tpc1; //Flash finder manager ::lightana::FlashFinderManager _mgr; ::lightana::FlashFinderManager _mgr_tpc0; @@ -217,7 +221,6 @@ private: double fPostWindow; double fPDFraction; - bool fDebug; art::ServiceHandle tfs; From 1d911cfdf273d688ab3be98de14232a0f615cb01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Tue, 7 Jul 2026 07:27:37 -0500 Subject: [PATCH 4/5] Add lightprop correction in default MC reco2 workflow --- .../JobConfigurations/standard/reco/config/workflow_reco2.fcl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/sbndcode/JobConfigurations/standard/reco/config/workflow_reco2.fcl b/sbndcode/JobConfigurations/standard/reco/config/workflow_reco2.fcl index fd2083d46..1654bf746 100755 --- a/sbndcode/JobConfigurations/standard/reco/config/workflow_reco2.fcl +++ b/sbndcode/JobConfigurations/standard/reco/config/workflow_reco2.fcl @@ -21,6 +21,7 @@ #include "SBNDCVNMapper.fcl" #include "blipreco_configs.fcl" #include "lightcalo_sbnd.fcl" +#include "sbnd_lightpropagationcorrection_config.fcl" BEGIN_PROLOG @@ -82,6 +83,7 @@ sbnd_reco2_producers:{ opt0finderSCE: @local::sbnd_opt0_finder_one_to_many tpcpmtbarycentermatching: @local::TPCPMTBarycenterMatchProducerMC tpcpmtbarycentermatchingSCE: @local::TPCPMTBarycenterMatchProducerSCEMC + lightpropagationcorrectionSCE: @local::LightPropagationCorrectionSCE ### Uncalibrated calorimetry producer for calibration caloskimmer caloskimCalorimetry: @local::caloskim_calorimetry @@ -136,6 +138,7 @@ sbnd_reco2_producer_sequence: [ #, fmatchoparaSCE , opt0finderSCE , tpcpmtbarycentermatchingSCE + , lightpropagationcorrectionSCE , cvn , cvnSCE , blipreco From 056ee388217cfa1534b53c24adeb49250ed011d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Tue, 21 Jul 2026 09:48:54 -0500 Subject: [PATCH 5/5] fcl refactoring --- sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd.fcl | 5 ++++- .../standard/caf/cafmakerjob_sbnd_data_sce.fcl | 1 - 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd.fcl b/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd.fcl index aa5aeafc7..ab521a814 100644 --- a/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd.fcl +++ b/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd.fcl @@ -156,9 +156,12 @@ physics.producers.cafmaker.FlashMatchSCECryoSuffixes: [""] physics.producers.cafmaker.TrackHitFillRRStartCut: 1000 physics.producers.cafmaker.TrackHitFillRREndCut: 1000 -# Include 3D barycenter flahs matching +# Include 3D barycenter flash matching physics.producers.cafmaker.TPCPMTBarycenterMatchLabel: "tpcpmtbarycentermatching" +# Include propagation correction on the flash timing +physics.producers.cafmaker.CorrectedOpFlashLabel: "lightpropagationcorrectionSCE" + # Include MC version of PMT Software Trigger physics.producers.cafmaker.SBNDSoftwareTriggerLabel: "pmtmetricmc" diff --git a/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd_data_sce.fcl b/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd_data_sce.fcl index 6aad2cdf6..552c9b465 100644 --- a/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd_data_sce.fcl +++ b/sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd_data_sce.fcl @@ -14,7 +14,6 @@ physics.producers.cafmaker.SBNDCRTTrackMatchLabel: "crttrackmatchingSCE" physics.producers.cafmaker.OpT0Label: "opt0finderSCE" physics.producers.cafmaker.TPCPMTBarycenterMatchLabel: "tpcpmtbarycentermatchingSCE" physics.producers.cafmaker.CVNLabel: "cvnSCE" -physics.producers.cafmaker.CorrectedOpFlashLabel: "lightpropagationcorrectionSCE" physics.producers.cafmaker.LightCaloLabel: "lightcaloSCE" physics.producers.cnnid.WireLabel: "sptpc2d:dnnsp"