Skip to content

Feature/acastill lightpropagation update #946

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
asanchezcastillo wants to merge 5 commits into
base: develop
Choose a base branch
from
Open

Feature/acastill lightpropagation update #946

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9c96580
Modify light propagation correction to include particle propagation time
asanchezcastillo Jul 6, 2026
c806e70
Add pdfraction as parameter
asanchezcastillo Jul 6, 2026
3e702ac
Clean code adding correctopflash function
asanchezcastillo Jul 7, 2026
1d911cf
Add lightprop correction in default MC reco2 workflow
asanchezcastillo Jul 7, 2026
056ee38
fcl refactoring
asanchezcastillo Jul 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down
1 change: 0 additions & 1 deletion sbndcode/JobConfigurations/standard/caf/cafmakerjob_sbnd_data_sce.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
3 changes: 3 additions & 0 deletions sbndcode/JobConfigurations/standard/reco/config/workflow_reco2.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
#include "SBNDCVNMapper.fcl"
#include "blipreco_configs.fcl"
#include "lightcalo_sbnd.fcl"
#include "sbnd_lightpropagationcorrection_config.fcl"

BEGIN_PROLOG

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -136,6 +138,7 @@ sbnd_reco2_producer_sequence: [
#, fmatchoparaSCE
, opt0finderSCE
, tpcpmtbarycentermatchingSCE
, lightpropagationcorrectionSCE
, cvn
, cvnSCE
, blipreco
Expand Down
10 changes: 8 additions & 2 deletions sbndcode/LightPropagationCorrection/LightPropagationCorrectionAna_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -95,7 +97,6 @@ void LightPropagationCorrectionAna::analyze(art::Event const& e)
return;
}


std::vector<art::Ptr<sbn::CorrectedOpFlashTiming>> correctedopflash_v;
art::fill_ptr_vector(correctedopflash_v, correctedopflash_h);

Expand All @@ -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();
}
}
Expand All @@ -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)
219 changes: 155 additions & 64 deletions sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,13 +22,19 @@ sbnd::LightPropagationCorrection::LightPropagationCorrection(fhicl::ParameterSet
fVGroupVUV( p.get<double>("VGroupVUV") ),
fNuScoreThreshold( p.get<double>("NuScoreThreshold") ),
fFMScoreThreshold( p.get<double>("FMScoreThreshold") ),
fMinHitPE ( p.get<int>("MinHitPE") ),
fPreWindow ( p.get<double>("PreWindow") ),
fPostWindow ( p.get<double>("PostWindow") ),
fPDFraction ( p.get<double>("PDFraction") ),
fDebug( p.get<bool>("Debug", false) )
//
// More initializers here.
{
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; opch<fNOpChannels; opch++){
auto pdCenter = fWireReadout.OpDetGeoFromOpChannel(opch).GetCenter();
Expand Down Expand Up @@ -113,8 +119,8 @@ void sbnd::LightPropagationCorrection::produce(art::Event & e)
//PFP to space points
art::FindManyP<recob::SpacePoint> pfp_sp_assns(pfpHandle, e, fSpacePointLabel);
//OpFlash to OpHit
art::FindManyP<recob::OpHit> flashToOpHitAssns_tpc0(opflashListHandle_tpc0, e, fOpFlashLabel_tpc0);
art::FindManyP<recob::OpHit> flashToOpHitAssns_tpc1(opflashListHandle_tpc1, e, fOpFlashLabel_tpc1);
flashToOpHitAssns_tpc0 = std::make_unique<art::FindManyP<recob::OpHit>>( opflashListHandle_tpc0, e, fOpFlashLabel_tpc0);

@tjones79 tjones79 Aug 11, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason these have to be declared as such?

I naively think these should be kept consistent with the above associations in how they are declared.

Could you move these back to local variables? You can pass them by const & into your new CorrectOpFlash() method. (Note: This is already done for slice_opt0finder_assns, slice_hit_assns etc... just above this line).

You would need to remove the corresponding declerations in the header file here

flashToOpHitAssns_tpc1 = std::make_unique<art::FindManyP<recob::OpHit>>(opflashListHandle_tpc1, e, fOpFlashLabel_tpc1);
// PFP Metadata
art::FindManyP<larpandoraobj::PFParticleMetadata> pfp_to_metadata(pfpHandle, e, fReco2Label);

Expand Down Expand Up @@ -231,63 +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<art::Ptr<recob::OpHit>> 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<recob::OpHit> newOpHitList;
std::vector<recob::OpHit> 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;
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));
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;
sbn::CorrectedOpFlashTiming correctedOpFlashTiming;
correctedOpFlashTiming.OpFlashT0 = originalFlashTime;
correctedOpFlashTiming.NuToFLight = (Zcenter/fSpeedOfLight)/1000;
correctedOpFlashTiming.NuToFCharge = (fRecoVz/fSpeedOfLight)/1000;
correctedOpFlashTiming.OpFlashT0Corrected = newFlashTime;
correctedOpFlashTimes->emplace_back(std::move(correctedOpFlashTiming));
}

sbn::CorrectedOpFlashTiming correctedOpFlashTiming;
CorrectOpFlash(flashFM[0], correctedOpFlashTiming);
correctedOpFlashTimes->emplace_back(std::move(correctedOpFlashTiming));
art::Ptr<sbn::CorrectedOpFlashTiming> 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));
Expand Down Expand Up @@ -396,6 +353,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);
Expand All @@ -408,13 +367,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<fSpacePointX.size(); sp++)
{
bool isInSameTPC = (fSpacePointX[sp] * _opDetX) > 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;
Expand All @@ -425,18 +388,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")
Comment on lines -428 to +393

@tjones79 tjones79 Aug 11, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A minor point, I think, but it looks like the local variable assignment const std::string pdType = fPDSMap.pdType(opdet); was dropped here, resulting in calling the function multiple times in the if/else if.

If restore the local variable then that can be avoided by simply doing: if(pdType=="pmt_coated" || pdType=="xarapuca_vuv") as in the original code

lightPropTime = lightPropTimeVIS;
else
throw std::runtime_error("LightPropagationCorrection: unexpected pdType '" + pdType + "' for opdet " + std::to_string(opdet));
Comment on lines -433 to -434

@tjones79 tjones79 Aug 11, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should restore the throw std::runtime_error to force the code to fail loudly if it sees something unexpected!

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;
}
}

Expand Down Expand Up @@ -515,12 +481,137 @@ 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) {
flash_hits_v.emplace_back(std::move(lite_hits_v.at(hitidx)));
}

return flash_hits_v;
}
}

double sbnd::LightPropagationCorrection::GetAverageParticlePropagationTime()
{
double sum = 0.0;
int n = 0;
for (size_t i=0; i<fSelectedChannelList.size(); i++) {
sum+= fParticlePropagationTimePerChannel[fSelectedChannelList[i]];
n++;
}
double average_prop_time = n ? sum / n : 0;
return average_prop_time;
return n ? sum / n : 0.0;

@tjones79 tjones79 Aug 11, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This return is unreachable and looks like a copy-and-paste error. Please remove!

}


double sbnd::LightPropagationCorrection::GetAveragePhotonPropagationTime()
{
double sum = 0.0;
int n = 0;
for (size_t i=0; i<fSelectedChannelList.size(); i++) {
sum+= fPhotonPropagationTimePerChannel[fSelectedChannelList[i]];
n++;
}
double average_prop_time = n ? sum / n : 0;
return average_prop_time;
return n ? sum / n : 0.0;

@tjones79 tjones79 Aug 11, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See previous comment

}

void sbnd::LightPropagationCorrection::GetSelectedChannelsFlash(
double flash_time,
::lightana::LiteOpHitArray_t ophit_list)
{
std::vector< std::tuple<double, double, size_t> > selected_hits;
double pe_sum = 0.0;
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]);

// Save channel ID
fSelectedChannelList.push_back(channel);

pe_count += pe;
if (pe_count / pe_sum > fPDFraction)
break;
}
}

void sbnd::LightPropagationCorrection::CorrectOpFlash(art::Ptr<recob::OpFlash> const& flash, sbn::CorrectedOpFlashTiming &correctedOpFlashTiming)
{
// Get the ophits associated to the flash
std::vector<art::Ptr<recob::OpHit>> 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<recob::OpHit> newOpHitList;
std::vector<recob::OpHit> 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);
}
}
Loading