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" 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 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..e064f715f 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc @@ -22,6 +22,10 @@ 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") ), + fPDFraction ( p.get("PDFraction") ), fDebug( p.get("Debug", false) ) // // More initializers here. @@ -29,6 +33,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; opch 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); @@ -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> 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; - 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 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)); @@ -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); @@ -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 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 +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") 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 +481,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 +489,129 @@ ::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; + 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 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 d0ae118fd..55b5944ff 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh @@ -142,10 +142,16 @@ 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(); + 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; @@ -176,6 +182,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 +215,12 @@ private: double fNuScoreThreshold; double fFMScoreThreshold; + // Parameters for the selected channel computation + int fMinHitPE; + double fPreWindow; + double fPostWindow; + double fPDFraction; + 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..ecfd6e2e8 100644 --- a/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl +++ b/sbndcode/LightPropagationCorrection/job/sbnd_lightpropagationcorrection_config.fcl @@ -35,6 +35,11 @@ LightPropagationCorrection: NuScoreThreshold: 0.1 FMScoreThreshold: 0.001 + MinHitPE: @local::FlashT0SelectedChannels.MinHitPE + PreWindow: @local::FlashT0SelectedChannels.PreWindow + PostWindow: @local::FlashT0SelectedChannels.PostWindow + PDFraction: @local::FlashT0SelectedChannels.PDFraction + FlashMatchingTool: "BarycenterFM" }