Feature/acastill lightpropagation update

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"
 

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"

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

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<art::Ptr<sbn::CorrectedOpFlashTiming>> 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)

sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc

@@ -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();
@@ -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);
+    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);
 
@@ -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));
@@ -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<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;
@@ -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,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;
+}
+
+
+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;
+}
+
+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);
+    }
+}