From 6063b9c05e7113d38712b9e0245862109bafa213 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Mon, 14 Jul 2025 04:06:16 -0500 Subject: [PATCH 1/8] Improve BFM with charge/light directionality --- .../TPCPMTBarycenterMatching/CMakeLists.txt | 4 + .../TPCPMTBarycenterMatching_module.cc | 511 +++++++++++++++--- ...sbnd_tpcpmt3dbarycentermatching_config.fcl | 37 +- 3 files changed, 478 insertions(+), 74 deletions(-) diff --git a/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt b/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt index 7158397dd..9a409d8f0 100644 --- a/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt +++ b/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt @@ -1,5 +1,6 @@ set( MODULE_LIBRARIES sbnobj::Common_Trigger + sbndcode_OpDetSim larcorealg::Geometry larcore::Geometry_Geometry_service lardataobj::RecoBase @@ -7,6 +8,8 @@ set( MODULE_LIBRARIES larcoreobj::SummaryData larana::OpticalDetector_OpHitFinder larreco::Calorimetry + larsim::PhotonPropagation + larsim::PhotonPropagation_PhotonVisibilityService_service lardata::Utilities nusimdata::SimulationBase nurandom::RandomUtils_NuRandomService_service @@ -31,6 +34,7 @@ set( MODULE_LIBRARIES ROOT::FFTW ROOT::Core ROOT::Tree + Eigen3::Eigen ) cet_build_plugin(TPCPMTBarycenterMatchProducer art::module SOURCE TPCPMTBarycenterMatching_module.cc LIBRARIES ${MODULE_LIBRARIES}) diff --git a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc index 8da5761d5..ceef2e434 100644 --- a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc +++ b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc @@ -39,7 +39,10 @@ #include "larcore/Geometry/Geometry.h" #include "larcore/Geometry/WireReadout.h" #include "larcorealg/Geometry/GeometryCore.h" +#include "larsim/PhotonPropagation/SemiAnalyticalModel.h" #include "larcore/CoreUtils/ServiceUtil.h" +#include "sbndcode/OpDetSim/sbndPDMapAlg.hh" + //Data product includes #include "lardataobj/RecoBase/OpHit.h" @@ -48,13 +51,19 @@ #include "lardataobj/RecoBase/Hit.h" #include "lardataobj/RecoBase/SpacePoint.h" #include "lardataobj/RecoBase/PFParticle.h" +#include "lardataobj/RecoBase/Cluster.h" #include "lardataobj/AnalysisBase/T0.h" #include "lardataobj/RawData/TriggerData.h" +#include "lardataobj/RecoBase/PFParticleMetadata.h" #include "sbnobj/Common/Reco/TPCPMTBarycenterMatch.h" //ROOT includes +#include +#include +#include +#include #include "TTree.h" -#include "TVector3.h" + #include // std::hypot(), std::abs(), std::sqrt() #include @@ -207,16 +216,34 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { void updateChargeVars(double sumCharge, TVector3 const& sumPos, TVector3 const& sumPosSqr, std::array const& triggerFlashCenter); ///< Update slice-level data members with charge and trigger match info void updateFlashVars(art::Ptr flash); ///< Update slice-level data members with best match info void updateMatchInfo(sbn::TPCPMTBarycenterMatch& matchInfo); ///< Update match product with slice-level data members - + void GetPCA(std::vector const& x, std::vector const& y, std::vector const& weight, std::vector& PCA ); ///< Get the PCA + double GetSliceCharge(const std::vector> &tpcHitsVec, const detinfo::DetectorPropertiesData det_prop, int tpc); + double GetFlashLight(double flash_pe, std::vector& total_dir_visibility, std::vector& total_ref_visibility); + void CreateOpHitList( std::vector> ophitlist, std::vector& ophit_z, std::vector& ophit_y, std::vector& ophit_weight); + // Input parameters std::vector fInputTags; ///< Suffix added onto fOpFlashLabel and fPandoraLabel, used by ICARUS for separate cryostat labels but could be empty std::vector fOpFlashesModuleLabel; ///< Label for PMT reconstruction products std::string fPandoraLabel; ///< Label for Pandora output products + std::string fTPCClusterLabel; bool fCollectionOnly; ///< Only use TPC spacepoints from the collection plane + double fDistanceCandidateFlashes; ///< Maximum distance between candidate flashes to be considered for matching (cm) + std::vector fCalAreaConst; /// Calibration area constants for wire plane + std::vector fSkipChannelList; + double fOpDetVUVEff; // Efficiencies for PMT detection + double fOpDetVISEff; // Efficiencies for PMT detection bool fVerbose; ///< Print extra info bool fFillMatchTree; ///< Fill an output TTree in the supplemental file bool fDo3DMatching; ///< Wether to perform the matching in 3D or 2D - + double fDistanceWeight; + double fAngleWeight; + std::vector fLightChargeRatioBounds; ///< Vector to store the distance between the barycenter of the charge and the barycenter of the light for each slice + double fXError; + double fYError; + double fZError; + double fAngleError; + std::vector fFlashVetoWindow; + int fDebugEvent; // Event-level data members int fRun; ///< Number of the run being processed int fEvent; ///< Number of the event being processed @@ -232,23 +259,41 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { double fChargeWidthY; ///< Weighted standard deviation of Y position of spacepoints (cm) double fChargeWidthZ; ///< Weighted standard deviation of Z position of spacepoints (cm) double fFlashTime; ///< Matched OpFlash time (us) + double fFlashTimeDebug; ///< Matched OpFlash time error (us) double fFlashPEs; ///< Brightness of matched flash (photoelectrons) double fFlashCenterX; ///< Flash position X obtained through eta_pmt curve (cm) double fFlashCenterY; ///< Weighted mean Y postion of hit PMTs (cm) double fFlashCenterZ; ///< Weighted mean Z postion of hit PMTs (cm) double fFlashWidthY; ///< Weighted standard deviation of Y postion of hit PMTs (cm) double fFlashWidthZ; ///< Weighted standard deviation of Z postion of hit PMTs (cm) + double fFlashLight; double fDeltaT; ///< | Matched flash time - charge T0 | when available (us) - double fDeltaX; ///< | Matched flash X center - charge X center | (cm) double fDeltaY; ///< | Matched flash Y center - charge Y center | (cm) double fDeltaZ; ///< | Matched flash Z center - charge Z center | (cm) double fRadius; ///< Hypotenuse of DeltaY and DeltaZ *parameter minimized by matching* (cm) - double fDeltaX_Trigger; ///< | Triggering flash X center - charge X center | (cm) double fDeltaY_Trigger; ///< | Triggering flash Y center - charge Y center | (cm) double fDeltaZ_Trigger; ///< | Triggering flash Z center - charge Z center | (cm) double fRadius_Trigger; ///< Hypotenuse of DeltaY_Trigger and DeltaZ_Trigger (cm) - + double fNuScore; ///< NuScore of the slice, if available, otherwise -9999 + double fNumElectrons; + double fLightChargeCosine; ///< Cosine of the angle between the charge PCA and the light PCA TTree* fMatchTree; ///< Tree to store all match information + TTree* fDebugTree; ///< Tree to store debug information + // Geometry service + geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); + opdet::sbndPDMapAlg fPDSMap; + + //Vector for PMT position + std::vector fOpDetID; + std::vector fOpDetType; + std::vector fOpDetX; + std::vector fOpDetY; + std::vector fOpDetZ; + + // Semi-analytical model for VUV and VIS light + std::unique_ptr _semi_model; + fhicl::ParameterSet _vuv_params; + fhicl::ParameterSet _vis_params; }; @@ -258,10 +303,24 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet // More initializers here. fOpFlashesModuleLabel(p.get>("OpFlashesModuleLabel")), fPandoraLabel(p.get("PandoraLabel")), + fTPCClusterLabel(p.get("TPCClusterLabel")), fCollectionOnly(p.get("CollectionOnly", true)), + fDistanceCandidateFlashes(p.get("DistanceCandidateFlashes")), // cm + fCalAreaConst(p.get>("CalAreaConst")), + fOpDetVUVEff (p.get("OpDetVUVEff")), + fOpDetVISEff (p.get("OpDetVISEff")), fVerbose(p.get("Verbose", false)), fFillMatchTree(p.get("FillMatchTree", false)), - fDo3DMatching(p.get("Do3DMatching",true)) + fDo3DMatching(p.get("Do3DMatching",true)), + fDistanceWeight(p.get("DistanceWeight")), + fAngleWeight(p.get("AngleWeight")), + fLightChargeRatioBounds(p.get>("LightChargeRatioBounds")), + fXError(p.get("XError")), // cm + fYError(p.get("YError")), // cm + fZError(p.get("ZError")), // cm + fAngleError(p.get("AngleError")), // deg + fFlashVetoWindow(p.get>("FlashVetoWindow")), // us + fDebugEvent(p.get("DebugEvent", -1)) // Event to debug, -1 means no debug event { // Call appropriate produces<>() functions here. @@ -274,9 +333,30 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet << "Module has been missconfigured. Number of OpFlashesModuleLabel must be 2!"; } - if ( fFillMatchTree ) { art::ServiceHandle tfs; - fMatchTree = tfs->make("matchTree","TPC Slice - OpFlash Matching Analysis"); + fDebugTree = tfs->make("debugTree","TPC Slice - OpFlash Matching Debug"); + //Event Info + fDebugTree->Branch("run", &fRun, "run/I" ); + fDebugTree->Branch("event", &fEvent, "event/I" ); + fDebugTree->Branch("cryo", &fTPC, "cryo/I" ); + fDebugTree->Branch("sliceNum", &fSliceNum, "sliceNum/I" ); + + fDebugTree->Branch("sliceX", &fChargeCenterX, "sliceX/d"); + fDebugTree->Branch("sliceY", &fChargeCenterY, "sliceY/d"); + fDebugTree->Branch("sliceZ", &fChargeCenterZ, "sliceZ/d"); + fDebugTree->Branch("sliceCharge", &fChargeTotal, "sliceCharge/d"); + fDebugTree->Branch("sliceElectrons", &fNumElectrons, "fNumElectrons/d" ); + fDebugTree->Branch("sliceNuScore", &fNuScore, "sliceNuScore/d"); + + fDebugTree->Branch("flashTime", &fFlashTimeDebug, "flashTime/d" ); + fDebugTree->Branch("flashPEs", &fFlashPEs, "flashPEs/d" ); + fDebugTree->Branch("flashLight", &fFlashLight, "flashLight/d" ); + fDebugTree->Branch("LightChargeCosine", &fLightChargeCosine, "LightChargeCosine/d" ); + + + if ( fFillMatchTree ) { + art::ServiceHandle tfs; + fMatchTree = tfs->make("matchTree","TPC Slice - OpFlash Matching Analysis"); //Event Info fMatchTree->Branch("run", &fRun, "run/I" ); @@ -293,6 +373,7 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet fMatchTree->Branch("chargeWidthX", &fChargeWidthX, "chargeWidthX/d" ); fMatchTree->Branch("chargeWidthY", &fChargeWidthY, "chargeWidthY/d" ); fMatchTree->Branch("chargeWidthZ", &fChargeWidthZ, "chargeWidthZ/d" ); + fMatchTree->Branch("nuScore", &fNuScore, "nuScore/d" ); //Matched Flash Info fMatchTree->Branch("flashTime", &fFlashTime, "flashTime/d" ); @@ -305,17 +386,33 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet //Match Quality Info fMatchTree->Branch("deltaT", &fDeltaT, "deltaT/d" ); - fMatchTree->Branch("deltaX", &fDeltaX, "deltaX/d" ); fMatchTree->Branch("deltaY", &fDeltaY, "deltaY/d" ); fMatchTree->Branch("deltaZ", &fDeltaZ, "deltaZ/d" ); fMatchTree->Branch("radius", &fRadius, "radius/d" ); - fMatchTree->Branch("deltaX_Trigger", &fDeltaX_Trigger, "deltaX_Trigger/d" ); fMatchTree->Branch("deltaZ_Trigger", &fDeltaZ_Trigger, "deltaZ_Trigger/d" ); fMatchTree->Branch("deltaY_Trigger", &fDeltaY_Trigger, "deltaY_Trigger/d" ); fMatchTree->Branch("radius_Trigger", &fRadius_Trigger, "dadius_Trigger/d" ); } //End MatchTree + //Fill the OpDet positions + for(unsigned int opch=0; opch("VUVHits"); + _vis_params = p.get("VIVHits"); + _semi_model = std::make_unique(_vuv_params, _vis_params, true, false); + } void TPCPMTBarycenterMatchProducer::produce(art::Event& e) @@ -324,68 +421,116 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) fEvent = e.id().event(); fRun = e.run(); + // Detector properties and clocks + auto const clock_data = art::ServiceHandle()->DataFor(e); + auto const det_prop = art::ServiceHandle()->DataFor(e, clock_data); + + //Initialize new data products auto matchInfoVector = std::make_unique< std::vector >(); art::PtrMaker< sbn::TPCPMTBarycenterMatch > const makeInfoPtr(e); auto sliceAssns = std::make_unique< art::Assns >(); auto flashAssns = std::make_unique< art::Assns >(); - /* ~~~~~~~~~~~~~~~~~~~~ Flash Section - * - * Here we gather the OpFlashes found in this cryostat and their OpHits - * We iterate through the flashes to identify a triggering flash - */ - std::array triggerFlashCenter ; - for ( size_t tpc=0; tpc<2 ; tpc ++) { - //Fetch the flashes and their associated hits, pointer vector needed for generating associations - art::Handle flashHandle = e.getHandle>(fOpFlashesModuleLabel[tpc]); - int nFlashes = (*flashHandle).size(); - triggerFlashCenter = {-9999., -9999., -9999.}; - double minTime = 99999.; - for (int i = 0; i < nFlashes; ++i) { - const recob::OpFlash &flash = (*flashHandle)[i]; - double _flashtime = flash.AbsTime(); - if ( std::abs(_flashtime)< minTime) { - minTime = std::abs(_flashtime); - triggerFlashCenter = {flash.XCenter(), flash.YCenter(), flash.ZCenter()}; - } + /* ~~~~~~~~~~~~~~~~~~~~ Flash Section + * + * Here we gather the OpFlashes found in this cryostat and their OpHits + * We iterate through the flashes to identify a triggering flash + */ + std::array triggerFlashCenter; + for ( size_t tpc=0; tpc<2 ; tpc ++) { + //Fetch the flashes and their associated hits, pointer vector needed for generating associations + art::Handle flashHandle = e.getHandle>(fOpFlashesModuleLabel[tpc]); + int nFlashes = (*flashHandle).size(); + triggerFlashCenter = {-9999., -9999., -9999.}; + double minTime = 99999.; + for (int i = 0; i < nFlashes; ++i) { + const recob::OpFlash &flash = (*flashHandle)[i]; + double _flashtime = flash.AbsTime(); + if ( std::abs(_flashtime)< minTime) { + minTime = std::abs(_flashtime); + triggerFlashCenter = {flash.XCenter(), flash.YCenter(), flash.ZCenter()}; + } + } } - } + // TODO (acastill): evaluate if there are two flashes in time coincidence and treat them as one. -/* ~~~~~~~~~~~~~~~~~~~~ TPC Section - * Here we start by gathering the Slices in the event - * For each slice, the charge centroid is first calculated - * Then we iterate through flashes to identify the best match flash - * If a triggering flash was found in this cyrostat, the barycenter distance to the triggering flash is also stored - */ + /* ~~~~~~~~~~~~~~~~~~~~ TPC Section + * Here we start by gathering the Slices in the event + * For each slice, the charge centroid is first calculated + * Then we iterate through flashes to identify the best match flash + * If a triggering flash was found in this cyrostat, the barycenter distance to the triggering flash is also stored + */ - //Fetch slices, TPC hits, and PFPs; pointer vector needed for generating associations - art::Handle const sliceHandle - = e.getHandle>(fPandoraLabel); + //Fetch slices, TPC hits, and PFPs; pointer vector needed for generating associations + ::art::Handle> sliceHandle; + e.getByLabel(fPandoraLabel, sliceHandle); art::FindManyP fmTPCHits(sliceHandle, e, fPandoraLabel); art::FindManyP fmPFPs(sliceHandle, e, fPandoraLabel); unsigned nSlices = (*sliceHandle).size(); + ::art::Handle> pfpHandle; + e.getByLabel(fPandoraLabel, pfpHandle); + //unsigned nPFPs = (*pfpHandle).size(); + art::FindManyP pfp_to_metadata(pfpHandle, e, fPandoraLabel); + art::FindManyP slice_pfp_assns (sliceHandle, e, fPandoraLabel); + + + // Hit to PFP assns + art::FindManyP fmHitsPFPs(pfpHandle, e, fPandoraLabel); + art::FindManyP fmClusterPfp(pfpHandle, e, fPandoraLabel); + ::art::Handle> clusterHandle; + e.getByLabel(fTPCClusterLabel, clusterHandle); + art::FindManyP cluster_hit_assns (clusterHandle, e, fPandoraLabel); + //For slice... for ( unsigned j = 0; j < nSlices; j++ ) { - //For each TPC - std::vector sliceMatchInfoVector; - std::vector> infoPtrVector; - std::vector> flashPtrVector; - fSliceNum = j; - const art::Ptr slicePtr { sliceHandle, j }; + //For each TPC + std::vector sliceMatchInfoVector; + std::vector> infoPtrVector; + std::vector> flashPtrVector; + fSliceNum = j; + const art::Ptr slicePtr { sliceHandle, j }; + // -----------------> This is to get the nuscore of the slice, only for debugging purposes <--------------------- // + //Slice to PFParticles association + fNuScore = -9999.; //Default value + //Vector for recob PFParticles + + std::vector> tpcHitsVec; + std::vector> pfpVect = slice_pfp_assns.at(j); + + // PFP Metadata + for(const art::Ptr &pfp : pfpVect){ + std::vector> cluster_v = fmClusterPfp.at(pfp.key()); + for(size_t i=0; i> hitVect = cluster_hit_assns.at(cluster_v[i].key()); + tpcHitsVec.insert(tpcHitsVec.end(), hitVect.begin(), hitVect.end()); + } + if(pfp->IsPrimary() && ( std::abs(pfp->PdgCode())==12 || std::abs(pfp->PdgCode())==14 ) ) + { + const std::vector> pfpMetaVec = pfp_to_metadata.at(pfp.key()); + for (auto const pfpMeta : pfpMetaVec) { + larpandoraobj::PFParticleMetadata::PropertiesMap propertiesMap = pfpMeta->GetPropertiesMap(); + fNuScore = propertiesMap.at("NuScore"); + } + } + } + for ( size_t tpc=0; tpc<2 ; tpc ++) { fTPC = tpc; InitializeSlice(); sbn::TPCPMTBarycenterMatch sliceMatchInfo; updateMatchInfo(sliceMatchInfo); - const std::vector> &tpcHitsVec = fmTPCHits.at(j); const std::vector> &pfpsVec = fmPFPs.at(j); art::FindOne f1SpacePoint(tpcHitsVec, e, fPandoraLabel); + std::vector hit_z; + std::vector hit_y; + std::vector hit_weight; + int nHits = tpcHitsVec.size(); int nPFPs = pfpsVec.size(); @@ -402,15 +547,34 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) TVector3 sumPos {0.,0.,0.}; TVector3 sumPosSqr {0.,0.,0.}; + size_t maxChargePlaneIdx=99999; + if(!fCollectionOnly) + { + std::vector PlaneCharge(3, 0.); // Vector to store the charge for each plane + //Loop to get the charge of each plane + for ( int k = 0; k < nHits; k++ ) { + if ( !f1SpacePoint.at(k).isValid() ) continue; + // If hit does not belong to correct tpc then skip + const art::Ptr &tpcHit = tpcHitsVec.at(k); + const recob::SpacePoint point = f1SpacePoint.at(k).ref(); + TVector3 const thisPoint = point.XYZ(); + if ((tpc == 0) == (thisPoint.X() > 0)) continue; // Skip if the point is not in the TPC we are considering + int plane = tpcHit->WireID().Plane; + PlaneCharge[plane] += tpcHit->Integral(); + } + // Idx to select the plane with the most charge + maxChargePlaneIdx = std::distance(PlaneCharge.begin(), std::max_element(PlaneCharge.begin(), PlaneCharge.end())); + } + //For hit... for ( int k = 0; k < nHits; k++ ) { // If hit does not belong to correct tpc then skip const art::Ptr &tpcHit = tpcHitsVec.at(k); //Only use hits with associated SpacePoints, and optionally only collection plane hits - if ( fCollectionOnly && tpcHit->SignalType() != geo::kCollection ) continue; if ( !f1SpacePoint.at(k).isValid() ) continue; - + if ( fCollectionOnly && tpcHit->WireID().Plane != 2 ) continue; + else if ( !fCollectionOnly && tpcHit->WireID().Plane != maxChargePlaneIdx) continue; // If not using collection plane only, skip hits not in the plane with the most charge const recob::SpacePoint point = f1SpacePoint.at(k).ref(); TVector3 const thisPoint = point.XYZ(); if ((tpc == 0) == (thisPoint.X() > 0)) continue; // Skip if the point is not in the TPC we are considering @@ -419,6 +583,10 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) sumCharge += thisCharge; sumPos += thisPoint * thisCharge; sumPosSqr += thisPointSqr * thisCharge; + // Store hit information for PCA + hit_y.push_back(thisPoint.Y()); + hit_z.push_back(thisPoint.Z()); + hit_weight.push_back(thisCharge); } //End for hit //No charge found in slice... @@ -431,27 +599,65 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) //Update charge variables updateChargeVars(sumCharge, sumPos, sumPosSqr, triggerFlashCenter); updateMatchInfo(sliceMatchInfo); + + std::vector ChargePCA = {0.,0.}; + GetPCA(hit_z, hit_y, hit_weight, ChargePCA); + + //Get the charge of the slice + double sliceCharge = this->GetSliceCharge(tpcHitsVec, det_prop, tpc); + fNumElectrons = sliceCharge; + //if(fNuScore>0.55) std::cout << " Event is " << fEvent << " Slice " << j << " in TPC " << tpc << " has barycenter at (" << fChargeCenterX << ", " << fChargeCenterY << ", " << fChargeCenterZ << ")" << " having used plane " << maxChargePlaneIdx << std::endl; + //Get the visibility map of the charge barycenter + geo::Point_t SliceXYZ = {fChargeCenterX, fChargeCenterY, fChargeCenterZ}; + std::cout << " Event " << fEvent<< " slice number " << j << " at tpc " << tpc << " with barycenter at (" << fChargeCenterX << ", " << fChargeCenterY << ", " << fChargeCenterZ << ")" << " with nuscore " << fNuScore << std::endl; + std::vector direct_visibility; + std::vector reflect_visibility; + _semi_model->detectedDirectVisibilities(direct_visibility, SliceXYZ); + _semi_model->detectedReflectedVisibilities(reflect_visibility, SliceXYZ); int matchIndex = -1; double minDistance = 1e6; double thisFlashCenterX, thisFlashCenterY, thisFlashCenterZ, thisDistance; - art::Handle flashHandle = e.getHandle>(fOpFlashesModuleLabel[tpc]); - int nFlashes = (*flashHandle).size(); + // --- Read Recob OpFlash + ::art::Handle> flashHandle; + e.getByLabel(fOpFlashesModuleLabel[tpc], flashHandle); + std::vector< art::Ptr > flashVect; + art::fill_ptr_vector(flashVect, flashHandle); + //OpHit OpFlash assns + art::FindManyP flash_ophit_assns(flashHandle, e, fOpFlashesModuleLabel[tpc]); + + //Vector to store the idxs of the candidate flashes + std::vector candidateFlashIdxs; + + int nFlashes = flashVect.size(); + + //Vector to store the distance between the barycenter of the charge and the barycenter of the light + std::vector chargeLightDistance; //For flash... for ( int m = 0; m < nFlashes; m++ ) { - const recob::OpFlash &flash = (*flashHandle).at(m); + auto & flash = flashVect[m]; //Find index of flash that minimizes barycenter distance in XYZ place - thisFlashCenterX = flash.XCenter(); - thisFlashCenterY = flash.YCenter(); - thisFlashCenterZ = flash.ZCenter(); + thisFlashCenterX = flash->XCenter(); + thisFlashCenterY = flash->YCenter(); + thisFlashCenterZ = flash->ZCenter(); + if(fDo3DMatching) thisDistance = std::hypot( (thisFlashCenterX - fChargeCenterX), (thisFlashCenterY - fChargeCenterY), (thisFlashCenterZ - fChargeCenterZ) ); else thisDistance = std::hypot( (thisFlashCenterY - fChargeCenterY), (thisFlashCenterZ - fChargeCenterZ) ); + chargeLightDistance.push_back(thisDistance); if ( thisDistance < minDistance ) { minDistance = thisDistance; matchIndex = m; } + // Fill the vector with the idxs of the flashes + if(thisDistanceTime() <<" has barycenter " << thisFlashCenterX << ", " << thisFlashCenterY << ", " << thisFlashCenterZ << " and distance " << thisDistance << " and PE " << flash->TotalPE() << std::endl; + } } //End for flash //No valid match found... @@ -461,6 +667,66 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) continue; } + // -----------------> Here we would end with the first pool of candidate flashes <----------------- + + double minChi2 = 1e6; + for(size_t i=0; iXCenter(); + double _currentFlashY = flash->YCenter(); + double _currentFlashZ = flash->ZCenter(); + //Get the ophits associated to this opflash + std::vector> ophitlist = flash_ophit_assns.at(flash.key()); + std::vector ophit_z; + std::vector ophit_y; + std::vector ophit_weight; + CreateOpHitList(ophitlist, ophit_z, ophit_y, ophit_weight); + + //std::cout << " Light PCA result for flash with idx " << idx << " with distance " << candidateFlashDistances[i] << std::endl; + double FlashLight = this->GetFlashLight(flash->TotalPE(), direct_visibility, reflect_visibility); + fFlashLight = FlashLight; + + fFlashTimeDebug = flash->Time(); + double lightChargeRatio = (FlashLight/sliceCharge); + std::vector LightPCA = {0.,0.}; + GetPCA(ophit_z, ophit_y, ophit_weight, LightPCA); + fLightChargeCosine = abs( ChargePCA[0]*LightPCA[0] + ChargePCA[1]*LightPCA[1] ) / ( std::hypot(ChargePCA[0], ChargePCA[1]) * std::hypot(LightPCA[0], LightPCA[1]) ); + double angle = std::acos(fLightChargeCosine)*(180/M_PI); + //double normalized_angle = angle/90.0; + //double normalized_distance = chargeLightDistance[idx]/fDistanceCandidateFlashes; + //double separation = fDistanceWeight*normalized_distance + fAngleWeight*normalized_angle; + double chi2; + if(fDo3DMatching) chi2 = std::pow(fChargeCenterX-_currentFlashX, 2)/std::pow(fXError, 2) + std::pow(fChargeCenterY-_currentFlashY, 2)/std::pow(fYError, 2) + std::pow(fChargeCenterZ-_currentFlashZ, 2)/std::pow(fZError, 2) + std::pow(angle, 2)/std::pow(fAngleError, 2); + else + chi2 = std::pow(fChargeCenterY-_currentFlashY, 2)/std::pow(fYError, 2) + std::pow(fChargeCenterZ-_currentFlashZ, 2)/std::pow(fZError, 2) + std::pow(angle, 2)/std::pow(fAngleError, 2); + + if(fNuScore>0.55) + { + std::cout << " flash " << i << " at time " << flash->Time() <<" has barycenter " << _currentFlashX << ", " << _currentFlashY << ", " << _currentFlashZ << " and distance " << chargeLightDistance[idx] << " with angle " << angle << " and charge light ratio " << lightChargeRatio << " and PE " << flash->TotalPE() <0.55) + { + std::cout << " --------> Candidate flash " << idx << " at time " << flash->Time() << " has LightChargeRatio " << lightChargeRatio <<" distance " << chargeLightDistance[idx] << " angle " << angle << " and chi2 " << chi2 << std::endl; + }*/ + if(lightChargeRatio < fLightChargeRatioBounds[0] || lightChargeRatio > fLightChargeRatioBounds[1]) continue; + if(chi2 < minChi2) + { + minChi2 = chi2; + matchIndex = idx; + } + fDebugTree->Fill(); + //double separation = fDistanceWeight*distance + fAngleWeight*angle; + } + + //std::cout << " Chosen flash is " << matchIndex << std::endl; + //Now we have a pool of flashes so we can have a look at the PCAs and the calorimetry + //Best match flash pointer unsigned unsignedMatchIndex = matchIndex; const art::Ptr flashPtr { flashHandle, unsignedMatchIndex }; @@ -476,7 +742,7 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) int maxFlashIdx=-1; int minNPes = -1000000000; - + for(int i=0; i(flashPtrVector.size()); i++) { int nPEs = flashPtrVector[i]->TotalPE(); @@ -486,12 +752,14 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) minNPes = nPEs; } } + if(maxFlashIdx!=-1) { sliceAssns->addSingle(infoPtrVector[maxFlashIdx], slicePtr); flashAssns->addSingle(infoPtrVector[maxFlashIdx], flashPtrVector[maxFlashIdx]); matchInfoVector->push_back(std::move(sliceMatchInfoVector[maxFlashIdx])); } + if(fNuScore>0.55) std::cout << " -----------------------> Selecting flash at time " << flashPtrVector[maxFlashIdx]->Time() << std::endl; } //End for slice //Store new products at the end of the event @@ -517,15 +785,15 @@ void TPCPMTBarycenterMatchProducer::InitializeSlice() { fFlashCenterZ = -9999.; fFlashWidthY = -9999.; fFlashWidthZ = -9999.; + fFlashLight = -9999.; fDeltaT = -9999.; - fDeltaX = -9999.; fDeltaY = -9999.; fDeltaZ = -9999.; fRadius = -9999.; - fDeltaX_Trigger = -9999.; fDeltaZ_Trigger = -9999.; fDeltaY_Trigger = -9999.; fRadius_Trigger = -9999.; + fNumElectrons = -9999.; } //End InitializeSlice() @@ -538,14 +806,8 @@ void TPCPMTBarycenterMatchProducer::updateChargeVars(double sumCharge, TVector3 fChargeWidthZ = std::sqrt( sumPosSqr[2]/sumCharge - (sumPos[2]/sumCharge)*(sumPos[2]/sumCharge) ); fChargeTotal = sumCharge; if ( triggerFlashCenter[1] != -9999. ) { - fDeltaX_Trigger = abs(triggerFlashCenter[0] - fChargeCenterX); - fDeltaY_Trigger = abs(triggerFlashCenter[1] - fChargeCenterY); - fDeltaZ_Trigger = abs(triggerFlashCenter[2] - fChargeCenterZ); - if(fDo3DMatching) - { - fRadius_Trigger = std::hypot(fDeltaX_Trigger, fDeltaY_Trigger, fDeltaZ_Trigger); - } - else + fDeltaY_Trigger = abs(triggerFlashCenter[0] - fChargeCenterY); + fDeltaZ_Trigger = abs(triggerFlashCenter[1] - fChargeCenterZ); fRadius_Trigger = std::hypot(fDeltaY_Trigger, fDeltaZ_Trigger); } } //End updateChargeVars() @@ -569,12 +831,6 @@ void TPCPMTBarycenterMatchProducer::updateFlashVars(art::Ptr fla if ( fChargeT0 != -9999 ) fDeltaT = abs(matchedTime - fChargeT0); fDeltaY = abs(matchedYCenter - fChargeCenterY); fDeltaZ = abs(matchedZCenter - fChargeCenterZ); - if( fDo3DMatching ) - { - fDeltaX = abs(matchedXCenter - fChargeCenterX); - fRadius = std::hypot(fDeltaX, fDeltaY, fDeltaZ); - } - else fRadius = std::hypot(fDeltaY, fDeltaZ); } //End updateFlashVars() @@ -594,5 +850,116 @@ void TPCPMTBarycenterMatchProducer::updateMatchInfo(sbn::TPCPMTBarycenterMatch& matchInfo.radius_Trigger = fRadius_Trigger; } //End updateMatchInfo() +void TPCPMTBarycenterMatchProducer::GetPCA(std::vector const&x, std::vector const&y, std::vector const&weight, std::vector & PCA ) { + + double weight_sum = std::accumulate(weight.begin(), weight.end(), 0.0); + //Create the matrix of points + size_t n = x.size(); + Eigen::MatrixXd points(n, 2); + for (size_t i = 0; i < n; ++i) { + points(i, 0) = x[i]; + points(i, 1) = y[i]; + } + //Get the weighted centroid + Eigen::Vector2d weighted_sum(0.0, 0.0); + for (size_t i = 0; i < n; ++i) { + weighted_sum(0) += weight[i] * x[i]; + weighted_sum(1) += weight[i] * y[i]; + } + Eigen::Vector2d centroid = weighted_sum / weight_sum; + //Center the points wrt centroid + Eigen::MatrixXd centered_points = points.rowwise() - centroid.transpose(); + //Get the weighted covariance matrix + Eigen::Matrix2d cov = Eigen::Matrix2d::Zero(); + for (size_t i = 0; i < n; ++i) { + Eigen::Vector2d pt = centered_points.row(i); + cov += std::abs(weight[i]) * (pt * pt.transpose()); + } + cov /= weight_sum; + // Obtener autovalores y autovectores + Eigen::SelfAdjointEigenSolver solver(cov); + Eigen::VectorXd eigvals = solver.eigenvalues(); + Eigen::Matrix2d eigvecs = solver.eigenvectors(); + + // Dirección principal (mayor autovalor) + int maxIndex; + eigvals.maxCoeff(&maxIndex); + Eigen::Vector2d principal_direction = eigvecs.col(maxIndex); + + PCA[0] = principal_direction(0); + PCA[1] = principal_direction(1); + +} + + + +double TPCPMTBarycenterMatchProducer::GetSliceCharge(const std::vector> &tpcHitsVec, const detinfo::DetectorPropertiesData det_prop, int tpc) { + + std::vector plane_charge{0.,0.,0.}; + std::vector plane_hits{0,0,0}; + for (size_t i=0; i < tpcHitsVec.size(); i++){ + auto hit = tpcHitsVec[i]; + if( hit->WireID().TPC != static_cast(tpc) ) continue; // Skip hits not in the TPC we are considering + auto drift_time = (hit->PeakTime() - 500)*0.5; // assuming TPC beam readout starts at 500 ticks, conversion = 0.5 us/tick + double atten_correction = std::exp(drift_time/det_prop.ElectronLifetime()); // exp(us/us) + auto hit_plane = hit->View(); + plane_charge.at(hit_plane) += hit->Integral()*atten_correction*(1/fCalAreaConst.at(hit_plane)); + plane_hits.at(hit_plane)++; + } + + //uint bestPlane = std::max_element(plane_charge.begin(), plane_charge.end()) - plane_charge.begin(); + uint bestHits = std::max_element(plane_hits.begin(), plane_hits.end()) - plane_hits.begin(); + + //double _mean_charge = (plane_charge[0] + plane_charge[1] + plane_charge[2])/3; + //double _max_charge = plane_charge.at(bestPlane); + double _comp_charge = plane_charge.at(bestHits); + //double _coll_charge = plane_charge[2]; + + return _comp_charge; +} + +void TPCPMTBarycenterMatchProducer::CreateOpHitList( std::vector> ophitlist, std::vector& ophit_z, std::vector& ophit_y, std::vector& ophit_weight) { + + std::map, double> footPrintMap; + for (size_t i=0; iOpChannel(); + int pos_z = static_cast(fOpDetZ[channel]); + int pos_y = static_cast(fOpDetY[channel]); + double weight = ophit->PE(); + footPrintMap[{pos_z, pos_y}] += weight; + } + + for (const auto& opdet : footPrintMap) { + int z = opdet.first.first; + int y = opdet.first.second; + int weight = opdet.second; + ophit_z.push_back(z); + ophit_y.push_back(y); + ophit_weight.push_back(weight*weight); + } + +} + +double TPCPMTBarycenterMatchProducer::GetFlashLight(double flash_pe, std::vector& dir_visibility, std::vector& ref_visibility) { + + double tot_visibility=0; + + for(size_t ch=0; ch Date: Mon, 14 Jul 2025 07:36:21 -0500 Subject: [PATCH 2/8] Update config fcl --- ...sbnd_tpcpmt3dbarycentermatching_config.fcl | 38 ++++--------------- 1 file changed, 7 insertions(+), 31 deletions(-) diff --git a/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl b/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl index cff283547..f50ea62b6 100644 --- a/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl +++ b/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl @@ -1,4 +1,6 @@ #include "opticalsimparameterisations_sbnd.fcl" +#include "digi_pmt_sbnd.fcl" +#include "calorimetry_sbnd.fcl" BEGIN_PROLOG @@ -6,7 +8,6 @@ TPCPMTBarycenterMatchProducer: { OpFlashesModuleLabel: ["opflashtpc0", "opflashtpc1"] PandoraLabel: "pandora" - TPCClusterLabel: "pandora" CollectionOnly: false Verbose: false FillMatchTree: false @@ -14,11 +15,9 @@ TPCPMTBarycenterMatchProducer: DistanceCandidateFlashes: 100 VUVHits: @local::sbnd_vuv_RS100cm_hits_parameterization VIVHits: @local::sbnd_vis_RS100cm_hits_parameterization - CalAreaConst: [ 0.0200906, 0.0200016, 0.0201293 ] - OpDetVUVEff: 0.032 - OpDetVISEff: 0.037 - DistanceWeight: 1 - AngleWeight: 0.45 + OpDetVUVEff: @local:sbnd_digipmt_alg.PMTCoatedVUVEff_tpc1 + OpDetVISEff: @local:sbnd_digipmt_alg.PMTUncoatedEff_tpc1 + CalAreaConst: @local:sbnd_calorimetryalgdata.CalAreaConstants LightChargeRatioBounds: [0.25, 3] XError: 12 YError: 23 @@ -28,31 +27,8 @@ TPCPMTBarycenterMatchProducer: module_type: "TPCPMTBarycenterMatchProducer" } -TPCPMTBarycenterMatchProducerSCE: -{ - OpFlashesModuleLabel: ["opflashtpc0", "opflashtpc1"] - PandoraLabel: "pandoraSCE" - TPCClusterLabel: "pandoraSCE" - CollectionOnly: false - Verbose: false - FillMatchTree: false - Do3DMatching: true - DistanceCandidateFlashes: 100 - VUVHits: @local::sbnd_vuv_RS100cm_hits_parameterization - VIVHits: @local::sbnd_vis_RS100cm_hits_parameterization - CalAreaConst: [ 0.0200906, 0.0200016, 0.0201293 ] - OpDetVUVEff: 0.032 - OpDetVISEff: 0.037 - DistanceWeight: 1 - AngleWeight: 0.45 - LightChargeRatioBounds: [0.25, 3] - XError: 12 - YError: 23 - ZError: 23 - AngleError: 27 - FlashVetoWindow: [-1500000,1500000] - module_type: "TPCPMTBarycenterMatchProducer" -} +TPCPMTBarycenterMatchProducerSCE: @local:TPCPMTBarycenterMatchProducer +TPCPMTBarycenterMatchProducerSCE.PandoraLabel: "pandoraSCE" END_PROLOG \ No newline at end of file From 34148277e86cbdc24e7bf0bcfa0afbdfbeb11a57 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Mon, 14 Jul 2025 07:36:46 -0500 Subject: [PATCH 3/8] Cleanup --- .../TPCPMTBarycenterMatching_module.cc | 52 ++----------------- 1 file changed, 4 insertions(+), 48 deletions(-) diff --git a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc index ceef2e434..39c0f6937 100644 --- a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc +++ b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc @@ -225,7 +225,6 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { std::vector fInputTags; ///< Suffix added onto fOpFlashLabel and fPandoraLabel, used by ICARUS for separate cryostat labels but could be empty std::vector fOpFlashesModuleLabel; ///< Label for PMT reconstruction products std::string fPandoraLabel; ///< Label for Pandora output products - std::string fTPCClusterLabel; bool fCollectionOnly; ///< Only use TPC spacepoints from the collection plane double fDistanceCandidateFlashes; ///< Maximum distance between candidate flashes to be considered for matching (cm) std::vector fCalAreaConst; /// Calibration area constants for wire plane @@ -235,8 +234,6 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { bool fVerbose; ///< Print extra info bool fFillMatchTree; ///< Fill an output TTree in the supplemental file bool fDo3DMatching; ///< Wether to perform the matching in 3D or 2D - double fDistanceWeight; - double fAngleWeight; std::vector fLightChargeRatioBounds; ///< Vector to store the distance between the barycenter of the charge and the barycenter of the light for each slice double fXError; double fYError; @@ -266,7 +263,6 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { double fFlashCenterZ; ///< Weighted mean Z postion of hit PMTs (cm) double fFlashWidthY; ///< Weighted standard deviation of Y postion of hit PMTs (cm) double fFlashWidthZ; ///< Weighted standard deviation of Z postion of hit PMTs (cm) - double fFlashLight; double fDeltaT; ///< | Matched flash time - charge T0 | when available (us) double fDeltaY; ///< | Matched flash Y center - charge Y center | (cm) double fDeltaZ; ///< | Matched flash Z center - charge Z center | (cm) @@ -275,10 +271,7 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { double fDeltaZ_Trigger; ///< | Triggering flash Z center - charge Z center | (cm) double fRadius_Trigger; ///< Hypotenuse of DeltaY_Trigger and DeltaZ_Trigger (cm) double fNuScore; ///< NuScore of the slice, if available, otherwise -9999 - double fNumElectrons; - double fLightChargeCosine; ///< Cosine of the angle between the charge PCA and the light PCA TTree* fMatchTree; ///< Tree to store all match information - TTree* fDebugTree; ///< Tree to store debug information // Geometry service geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); opdet::sbndPDMapAlg fPDSMap; @@ -303,7 +296,6 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet // More initializers here. fOpFlashesModuleLabel(p.get>("OpFlashesModuleLabel")), fPandoraLabel(p.get("PandoraLabel")), - fTPCClusterLabel(p.get("TPCClusterLabel")), fCollectionOnly(p.get("CollectionOnly", true)), fDistanceCandidateFlashes(p.get("DistanceCandidateFlashes")), // cm fCalAreaConst(p.get>("CalAreaConst")), @@ -312,8 +304,6 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet fVerbose(p.get("Verbose", false)), fFillMatchTree(p.get("FillMatchTree", false)), fDo3DMatching(p.get("Do3DMatching",true)), - fDistanceWeight(p.get("DistanceWeight")), - fAngleWeight(p.get("AngleWeight")), fLightChargeRatioBounds(p.get>("LightChargeRatioBounds")), fXError(p.get("XError")), // cm fYError(p.get("YError")), // cm @@ -333,27 +323,7 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet << "Module has been missconfigured. Number of OpFlashesModuleLabel must be 2!"; } - art::ServiceHandle tfs; - fDebugTree = tfs->make("debugTree","TPC Slice - OpFlash Matching Debug"); - //Event Info - fDebugTree->Branch("run", &fRun, "run/I" ); - fDebugTree->Branch("event", &fEvent, "event/I" ); - fDebugTree->Branch("cryo", &fTPC, "cryo/I" ); - fDebugTree->Branch("sliceNum", &fSliceNum, "sliceNum/I" ); - - fDebugTree->Branch("sliceX", &fChargeCenterX, "sliceX/d"); - fDebugTree->Branch("sliceY", &fChargeCenterY, "sliceY/d"); - fDebugTree->Branch("sliceZ", &fChargeCenterZ, "sliceZ/d"); - fDebugTree->Branch("sliceCharge", &fChargeTotal, "sliceCharge/d"); - fDebugTree->Branch("sliceElectrons", &fNumElectrons, "fNumElectrons/d" ); - fDebugTree->Branch("sliceNuScore", &fNuScore, "sliceNuScore/d"); - - fDebugTree->Branch("flashTime", &fFlashTimeDebug, "flashTime/d" ); - fDebugTree->Branch("flashPEs", &fFlashPEs, "flashPEs/d" ); - fDebugTree->Branch("flashLight", &fFlashLight, "flashLight/d" ); - fDebugTree->Branch("LightChargeCosine", &fLightChargeCosine, "LightChargeCosine/d" ); - - + art::ServiceHandle tfs; if ( fFillMatchTree ) { art::ServiceHandle tfs; fMatchTree = tfs->make("matchTree","TPC Slice - OpFlash Matching Analysis"); @@ -481,7 +451,7 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) art::FindManyP fmHitsPFPs(pfpHandle, e, fPandoraLabel); art::FindManyP fmClusterPfp(pfpHandle, e, fPandoraLabel); ::art::Handle> clusterHandle; - e.getByLabel(fTPCClusterLabel, clusterHandle); + e.getByLabel(fPandoraLabel, clusterHandle); art::FindManyP cluster_hit_assns (clusterHandle, e, fPandoraLabel); @@ -605,11 +575,8 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) //Get the charge of the slice double sliceCharge = this->GetSliceCharge(tpcHitsVec, det_prop, tpc); - fNumElectrons = sliceCharge; - //if(fNuScore>0.55) std::cout << " Event is " << fEvent << " Slice " << j << " in TPC " << tpc << " has barycenter at (" << fChargeCenterX << ", " << fChargeCenterY << ", " << fChargeCenterZ << ")" << " having used plane " << maxChargePlaneIdx << std::endl; //Get the visibility map of the charge barycenter geo::Point_t SliceXYZ = {fChargeCenterX, fChargeCenterY, fChargeCenterZ}; - std::cout << " Event " << fEvent<< " slice number " << j << " at tpc " << tpc << " with barycenter at (" << fChargeCenterX << ", " << fChargeCenterY << ", " << fChargeCenterZ << ")" << " with nuscore " << fNuScore << std::endl; std::vector direct_visibility; std::vector reflect_visibility; _semi_model->detectedDirectVisibilities(direct_visibility, SliceXYZ); @@ -673,7 +640,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) for(size_t i=0; iXCenter(); @@ -686,19 +652,13 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) std::vector ophit_weight; CreateOpHitList(ophitlist, ophit_z, ophit_y, ophit_weight); - //std::cout << " Light PCA result for flash with idx " << idx << " with distance " << candidateFlashDistances[i] << std::endl; double FlashLight = this->GetFlashLight(flash->TotalPE(), direct_visibility, reflect_visibility); - fFlashLight = FlashLight; - fFlashTimeDebug = flash->Time(); double lightChargeRatio = (FlashLight/sliceCharge); std::vector LightPCA = {0.,0.}; GetPCA(ophit_z, ophit_y, ophit_weight, LightPCA); - fLightChargeCosine = abs( ChargePCA[0]*LightPCA[0] + ChargePCA[1]*LightPCA[1] ) / ( std::hypot(ChargePCA[0], ChargePCA[1]) * std::hypot(LightPCA[0], LightPCA[1]) ); - double angle = std::acos(fLightChargeCosine)*(180/M_PI); - //double normalized_angle = angle/90.0; - //double normalized_distance = chargeLightDistance[idx]/fDistanceCandidateFlashes; - //double separation = fDistanceWeight*normalized_distance + fAngleWeight*normalized_angle; + double LightChargeCosine = abs( ChargePCA[0]*LightPCA[0] + ChargePCA[1]*LightPCA[1] ) / ( std::hypot(ChargePCA[0], ChargePCA[1]) * std::hypot(LightPCA[0], LightPCA[1]) ); + double angle = std::acos(LightChargeCosine)*(180/M_PI); double chi2; if(fDo3DMatching) chi2 = std::pow(fChargeCenterX-_currentFlashX, 2)/std::pow(fXError, 2) + std::pow(fChargeCenterY-_currentFlashY, 2)/std::pow(fYError, 2) + std::pow(fChargeCenterZ-_currentFlashZ, 2)/std::pow(fZError, 2) + std::pow(angle, 2)/std::pow(fAngleError, 2); else @@ -720,8 +680,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) minChi2 = chi2; matchIndex = idx; } - fDebugTree->Fill(); - //double separation = fDistanceWeight*distance + fAngleWeight*angle; } //std::cout << " Chosen flash is " << matchIndex << std::endl; @@ -785,7 +743,6 @@ void TPCPMTBarycenterMatchProducer::InitializeSlice() { fFlashCenterZ = -9999.; fFlashWidthY = -9999.; fFlashWidthZ = -9999.; - fFlashLight = -9999.; fDeltaT = -9999.; fDeltaY = -9999.; fDeltaZ = -9999.; @@ -793,7 +750,6 @@ void TPCPMTBarycenterMatchProducer::InitializeSlice() { fDeltaZ_Trigger = -9999.; fDeltaY_Trigger = -9999.; fRadius_Trigger = -9999.; - fNumElectrons = -9999.; } //End InitializeSlice() From b8074206ae318376e412f03f708868f041303cb3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Mon, 14 Jul 2025 08:16:25 -0500 Subject: [PATCH 4/8] Cleanup --- .../TPCPMTBarycenterMatching_module.cc | 51 ++++--------------- 1 file changed, 10 insertions(+), 41 deletions(-) diff --git a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc index 39c0f6937..2decc82dd 100644 --- a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc +++ b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc @@ -240,7 +240,6 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { double fZError; double fAngleError; std::vector fFlashVetoWindow; - int fDebugEvent; // Event-level data members int fRun; ///< Number of the run being processed int fEvent; ///< Number of the event being processed @@ -270,7 +269,6 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { double fDeltaY_Trigger; ///< | Triggering flash Y center - charge Y center | (cm) double fDeltaZ_Trigger; ///< | Triggering flash Z center - charge Z center | (cm) double fRadius_Trigger; ///< Hypotenuse of DeltaY_Trigger and DeltaZ_Trigger (cm) - double fNuScore; ///< NuScore of the slice, if available, otherwise -9999 TTree* fMatchTree; ///< Tree to store all match information // Geometry service geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); @@ -310,7 +308,6 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet fZError(p.get("ZError")), // cm fAngleError(p.get("AngleError")), // deg fFlashVetoWindow(p.get>("FlashVetoWindow")), // us - fDebugEvent(p.get("DebugEvent", -1)) // Event to debug, -1 means no debug event { // Call appropriate produces<>() functions here. @@ -343,7 +340,6 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet fMatchTree->Branch("chargeWidthX", &fChargeWidthX, "chargeWidthX/d" ); fMatchTree->Branch("chargeWidthY", &fChargeWidthY, "chargeWidthY/d" ); fMatchTree->Branch("chargeWidthZ", &fChargeWidthZ, "chargeWidthZ/d" ); - fMatchTree->Branch("nuScore", &fNuScore, "nuScore/d" ); //Matched Flash Info fMatchTree->Branch("flashTime", &fFlashTime, "flashTime/d" ); @@ -442,7 +438,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) unsigned nSlices = (*sliceHandle).size(); ::art::Handle> pfpHandle; e.getByLabel(fPandoraLabel, pfpHandle); - //unsigned nPFPs = (*pfpHandle).size(); art::FindManyP pfp_to_metadata(pfpHandle, e, fPandoraLabel); art::FindManyP slice_pfp_assns (sliceHandle, e, fPandoraLabel); @@ -465,27 +460,17 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) const art::Ptr slicePtr { sliceHandle, j }; // -----------------> This is to get the nuscore of the slice, only for debugging purposes <--------------------- // //Slice to PFParticles association - fNuScore = -9999.; //Default value //Vector for recob PFParticles - std::vector> tpcHitsVec; std::vector> pfpVect = slice_pfp_assns.at(j); - // PFP Metadata + // Get the hits associated to the PFParticles in this slice for(const art::Ptr &pfp : pfpVect){ std::vector> cluster_v = fmClusterPfp.at(pfp.key()); for(size_t i=0; i> hitVect = cluster_hit_assns.at(cluster_v[i].key()); tpcHitsVec.insert(tpcHitsVec.end(), hitVect.begin(), hitVect.end()); } - if(pfp->IsPrimary() && ( std::abs(pfp->PdgCode())==12 || std::abs(pfp->PdgCode())==14 ) ) - { - const std::vector> pfpMetaVec = pfp_to_metadata.at(pfp.key()); - for (auto const pfpMeta : pfpMetaVec) { - larpandoraobj::PFParticleMetadata::PropertiesMap propertiesMap = pfpMeta->GetPropertiesMap(); - fNuScore = propertiesMap.at("NuScore"); - } - } } for ( size_t tpc=0; tpc<2 ; tpc ++) { @@ -591,16 +576,14 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) e.getByLabel(fOpFlashesModuleLabel[tpc], flashHandle); std::vector< art::Ptr > flashVect; art::fill_ptr_vector(flashVect, flashHandle); + //OpHit OpFlash assns art::FindManyP flash_ophit_assns(flashHandle, e, fOpFlashesModuleLabel[tpc]); //Vector to store the idxs of the candidate flashes std::vector candidateFlashIdxs; - int nFlashes = flashVect.size(); - //Vector to store the distance between the barycenter of the charge and the barycenter of the light - std::vector chargeLightDistance; //For flash... for ( int m = 0; m < nFlashes; m++ ) { auto & flash = flashVect[m]; @@ -612,7 +595,7 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) if(fDo3DMatching) thisDistance = std::hypot( (thisFlashCenterX - fChargeCenterX), (thisFlashCenterY - fChargeCenterY), (thisFlashCenterZ - fChargeCenterZ) ); else thisDistance = std::hypot( (thisFlashCenterY - fChargeCenterY), (thisFlashCenterZ - fChargeCenterZ) ); - chargeLightDistance.push_back(thisDistance); + if ( thisDistance < minDistance ) { minDistance = thisDistance; matchIndex = m; @@ -621,10 +604,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) if(thisDistanceTime() <<" has barycenter " << thisFlashCenterX << ", " << thisFlashCenterY << ", " << thisFlashCenterZ << " and distance " << thisDistance << " and PE " << flash->TotalPE() << std::endl; - } } //End for flash //No valid match found... @@ -664,16 +643,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) else chi2 = std::pow(fChargeCenterY-_currentFlashY, 2)/std::pow(fYError, 2) + std::pow(fChargeCenterZ-_currentFlashZ, 2)/std::pow(fZError, 2) + std::pow(angle, 2)/std::pow(fAngleError, 2); - if(fNuScore>0.55) - { - std::cout << " flash " << i << " at time " << flash->Time() <<" has barycenter " << _currentFlashX << ", " << _currentFlashY << ", " << _currentFlashZ << " and distance " << chargeLightDistance[idx] << " with angle " << angle << " and charge light ratio " << lightChargeRatio << " and PE " << flash->TotalPE() <0.55) - { - std::cout << " --------> Candidate flash " << idx << " at time " << flash->Time() << " has LightChargeRatio " << lightChargeRatio <<" distance " << chargeLightDistance[idx] << " angle " << angle << " and chi2 " << chi2 << std::endl; - }*/ if(lightChargeRatio < fLightChargeRatioBounds[0] || lightChargeRatio > fLightChargeRatioBounds[1]) continue; if(chi2 < minChi2) { @@ -682,9 +651,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) } } - //std::cout << " Chosen flash is " << matchIndex << std::endl; - //Now we have a pool of flashes so we can have a look at the PCAs and the calorimetry - //Best match flash pointer unsigned unsignedMatchIndex = matchIndex; const art::Ptr flashPtr { flashHandle, unsignedMatchIndex }; @@ -717,7 +683,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) flashAssns->addSingle(infoPtrVector[maxFlashIdx], flashPtrVector[maxFlashIdx]); matchInfoVector->push_back(std::move(sliceMatchInfoVector[maxFlashIdx])); } - if(fNuScore>0.55) std::cout << " -----------------------> Selecting flash at time " << flashPtrVector[maxFlashIdx]->Time() << std::endl; } //End for slice //Store new products at the end of the event @@ -816,15 +781,18 @@ void TPCPMTBarycenterMatchProducer::GetPCA(std::vector const&x, std::vec points(i, 0) = x[i]; points(i, 1) = y[i]; } + //Get the weighted centroid Eigen::Vector2d weighted_sum(0.0, 0.0); for (size_t i = 0; i < n; ++i) { weighted_sum(0) += weight[i] * x[i]; weighted_sum(1) += weight[i] * y[i]; } - Eigen::Vector2d centroid = weighted_sum / weight_sum; + //Center the points wrt centroid + Eigen::Vector2d centroid = weighted_sum / weight_sum; Eigen::MatrixXd centered_points = points.rowwise() - centroid.transpose(); + //Get the weighted covariance matrix Eigen::Matrix2d cov = Eigen::Matrix2d::Zero(); for (size_t i = 0; i < n; ++i) { @@ -832,12 +800,13 @@ void TPCPMTBarycenterMatchProducer::GetPCA(std::vector const&x, std::vec cov += std::abs(weight[i]) * (pt * pt.transpose()); } cov /= weight_sum; - // Obtener autovalores y autovectores + + // Get eigenvalues / eigenvectors Eigen::SelfAdjointEigenSolver solver(cov); Eigen::VectorXd eigvals = solver.eigenvalues(); Eigen::Matrix2d eigvecs = solver.eigenvectors(); - // Dirección principal (mayor autovalor) + // Get Principal component (largest eigenvalue) int maxIndex; eigvals.maxCoeff(&maxIndex); Eigen::Vector2d principal_direction = eigvecs.col(maxIndex); From efb42abad134ea345d95442ff5a1062264bcee6b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Mon, 14 Jul 2025 08:29:07 -0500 Subject: [PATCH 5/8] Cleanup --- .../TPCPMTBarycenterMatching_module.cc | 54 ++++++++----------- 1 file changed, 23 insertions(+), 31 deletions(-) diff --git a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc index 2decc82dd..ec504327c 100644 --- a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc +++ b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc @@ -432,21 +432,22 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) //Fetch slices, TPC hits, and PFPs; pointer vector needed for generating associations ::art::Handle> sliceHandle; - e.getByLabel(fPandoraLabel, sliceHandle); - art::FindManyP fmTPCHits(sliceHandle, e, fPandoraLabel); - art::FindManyP fmPFPs(sliceHandle, e, fPandoraLabel); + e.getByLabel(fPandoraLabel, sliceHandle); + + // Slice to PFP assns unsigned nSlices = (*sliceHandle).size(); ::art::Handle> pfpHandle; e.getByLabel(fPandoraLabel, pfpHandle); - art::FindManyP pfp_to_metadata(pfpHandle, e, fPandoraLabel); - art::FindManyP slice_pfp_assns (sliceHandle, e, fPandoraLabel); + // PFP to metadata assns + art::FindManyP slice_pfp_assns(sliceHandle, e, fPandoraLabel); // Hit to PFP assns - art::FindManyP fmHitsPFPs(pfpHandle, e, fPandoraLabel); - art::FindManyP fmClusterPfp(pfpHandle, e, fPandoraLabel); + art::FindManyP pfp_cluster_assns(pfpHandle, e, fPandoraLabel); ::art::Handle> clusterHandle; e.getByLabel(fPandoraLabel, clusterHandle); + + // Cluster to Hit assns art::FindManyP cluster_hit_assns (clusterHandle, e, fPandoraLabel); @@ -458,44 +459,40 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) std::vector> flashPtrVector; fSliceNum = j; const art::Ptr slicePtr { sliceHandle, j }; - // -----------------> This is to get the nuscore of the slice, only for debugging purposes <--------------------- // - //Slice to PFParticles association //Vector for recob PFParticles std::vector> tpcHitsVec; std::vector> pfpVect = slice_pfp_assns.at(j); // Get the hits associated to the PFParticles in this slice for(const art::Ptr &pfp : pfpVect){ - std::vector> cluster_v = fmClusterPfp.at(pfp.key()); + std::vector> cluster_v = pfp_cluster_assns.at(pfp.key()); for(size_t i=0; i> hitVect = cluster_hit_assns.at(cluster_v[i].key()); tpcHitsVec.insert(tpcHitsVec.end(), hitVect.begin(), hitVect.end()); } } + int nPFPs = pfpVect.size(); + //Retrieve Pandora's T0 for this slice if available, same for every PFP in slice so we only need one + if ( nPFPs != 0 ) { + art::FindOne f1T0( {pfpVect.at(0)}, e, fPandoraLabel); + if ( f1T0.at(0).isValid() ) { + fChargeT0 = f1T0.at(0).ref().Time() / 1e3; + } + } + for ( size_t tpc=0; tpc<2 ; tpc ++) { fTPC = tpc; InitializeSlice(); sbn::TPCPMTBarycenterMatch sliceMatchInfo; updateMatchInfo(sliceMatchInfo); - - const std::vector> &pfpsVec = fmPFPs.at(j); art::FindOne f1SpacePoint(tpcHitsVec, e, fPandoraLabel); std::vector hit_z; std::vector hit_y; std::vector hit_weight; - - int nHits = tpcHitsVec.size(); - int nPFPs = pfpsVec.size(); - //Retrieve Pandora's T0 for this slice if available, same for every PFP in slice so we only need one - if ( nPFPs != 0 ) { - art::FindOne f1T0( {pfpsVec.at(0)}, e, fPandoraLabel); - if ( f1T0.at(0).isValid() ) { - fChargeT0 = f1T0.at(0).ref().Time() / 1e3; - } - } + int nHits = tpcHitsVec.size(); double thisCharge; double sumCharge = 0.; @@ -503,8 +500,9 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) TVector3 sumPosSqr {0.,0.,0.}; size_t maxChargePlaneIdx=99999; - if(!fCollectionOnly) - { + + // If we are not using collection plane only, we need to find the plane with the most charge + if(!fCollectionOnly){ std::vector PlaneCharge(3, 0.); // Vector to store the charge for each plane //Loop to get the charge of each plane for ( int k = 0; k < nHits; k++ ) { @@ -832,13 +830,8 @@ double TPCPMTBarycenterMatchProducer::GetSliceCharge(const std::vector Date: Mon, 14 Jul 2025 08:44:57 -0500 Subject: [PATCH 6/8] Fix typo --- .../TPCPMTBarycenterMatching_module.cc | 10 +++++----- .../sbnd_tpcpmt3dbarycentermatching_config.fcl | 18 +++++++++--------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc index ec504327c..e38f5a720 100644 --- a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc +++ b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc @@ -299,15 +299,15 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet fCalAreaConst(p.get>("CalAreaConst")), fOpDetVUVEff (p.get("OpDetVUVEff")), fOpDetVISEff (p.get("OpDetVISEff")), - fVerbose(p.get("Verbose", false)), - fFillMatchTree(p.get("FillMatchTree", false)), - fDo3DMatching(p.get("Do3DMatching",true)), + fVerbose(p.get("Verbose")), + fFillMatchTree(p.get("FillMatchTree")), + fDo3DMatching(p.get("Do3DMatching")), fLightChargeRatioBounds(p.get>("LightChargeRatioBounds")), fXError(p.get("XError")), // cm fYError(p.get("YError")), // cm fZError(p.get("ZError")), // cm fAngleError(p.get("AngleError")), // deg - fFlashVetoWindow(p.get>("FlashVetoWindow")), // us + fFlashVetoWindow(p.get>("FlashVetoWindow")) // us { // Call appropriate produces<>() functions here. @@ -328,7 +328,7 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet //Event Info fMatchTree->Branch("run", &fRun, "run/I" ); fMatchTree->Branch("event", &fEvent, "event/I" ); - fMatchTree->Branch("cryo", &fTPC, "cryo/I" ); + fMatchTree->Branch("cryo", &fTPC, "cryo/I" ); fMatchTree->Branch("sliceNum", &fSliceNum, "sliceNum/I" ); //Charge Info diff --git a/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl b/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl index f50ea62b6..d6f622c8a 100644 --- a/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl +++ b/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl @@ -15,19 +15,19 @@ TPCPMTBarycenterMatchProducer: DistanceCandidateFlashes: 100 VUVHits: @local::sbnd_vuv_RS100cm_hits_parameterization VIVHits: @local::sbnd_vis_RS100cm_hits_parameterization - OpDetVUVEff: @local:sbnd_digipmt_alg.PMTCoatedVUVEff_tpc1 - OpDetVISEff: @local:sbnd_digipmt_alg.PMTUncoatedEff_tpc1 - CalAreaConst: @local:sbnd_calorimetryalgdata.CalAreaConstants + OpDetVUVEff: @local::sbnd_digipmt_alg.PMTCoatedVUVEff_tpc1 + OpDetVISEff: @local::sbnd_digipmt_alg.PMTUncoatedEff_tpc1 + CalAreaConst: @local::sbnd_calorimetryalgdata.CalAreaConstants LightChargeRatioBounds: [0.25, 3] - XError: 12 - YError: 23 - ZError: 23 - AngleError: 27 - FlashVetoWindow: [-1500000,1500000] + XError: 12 // cm + YError: 23 // cm + ZError: 23 // cm + AngleError: 27 // deg + FlashVetoWindow: [-1500000,1500000] // ns module_type: "TPCPMTBarycenterMatchProducer" } -TPCPMTBarycenterMatchProducerSCE: @local:TPCPMTBarycenterMatchProducer +TPCPMTBarycenterMatchProducerSCE: @local::TPCPMTBarycenterMatchProducer TPCPMTBarycenterMatchProducerSCE.PandoraLabel: "pandoraSCE" From eefb59e6ce2d08ffff788af288a344fb7062e8ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= Date: Tue, 2 Sep 2025 08:11:46 -0500 Subject: [PATCH 7/8] Add score to the tpcpmt bfm --- .../TPCPMTBarycenterMatching_module.cc | 24 ++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc index e38f5a720..a9f260fe3 100644 --- a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc +++ b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc @@ -225,6 +225,7 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { std::vector fInputTags; ///< Suffix added onto fOpFlashLabel and fPandoraLabel, used by ICARUS for separate cryostat labels but could be empty std::vector fOpFlashesModuleLabel; ///< Label for PMT reconstruction products std::string fPandoraLabel; ///< Label for Pandora output products + std::string fOpT0Label; bool fCollectionOnly; ///< Only use TPC spacepoints from the collection plane double fDistanceCandidateFlashes; ///< Maximum distance between candidate flashes to be considered for matching (cm) std::vector fCalAreaConst; /// Calibration area constants for wire plane @@ -266,6 +267,9 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { double fDeltaY; ///< | Matched flash Y center - charge Y center | (cm) double fDeltaZ; ///< | Matched flash Z center - charge Z center | (cm) double fRadius; ///< Hypotenuse of DeltaY and DeltaZ *parameter minimized by matching* (cm) + double fChi2; ///< Chi2 to be minimized when matching flash to slice (dimensionless) + double fScore; ///< Score to be maximized when matching flash to slice (dimensionless) + double fAngle; ///< Angle between charge PCA and light PCA (degrees) double fDeltaY_Trigger; ///< | Triggering flash Y center - charge Y center | (cm) double fDeltaZ_Trigger; ///< | Triggering flash Z center - charge Z center | (cm) double fRadius_Trigger; ///< Hypotenuse of DeltaY_Trigger and DeltaZ_Trigger (cm) @@ -355,9 +359,12 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet fMatchTree->Branch("deltaY", &fDeltaY, "deltaY/d" ); fMatchTree->Branch("deltaZ", &fDeltaZ, "deltaZ/d" ); fMatchTree->Branch("radius", &fRadius, "radius/d" ); + fMatchTree->Branch("angle", &fAngle, "angle/d" ); + fMatchTree->Branch("chi2", &fChi2, "chi2/d" ); + fMatchTree->Branch("score", &fScore, "score/d" ); fMatchTree->Branch("deltaZ_Trigger", &fDeltaZ_Trigger, "deltaZ_Trigger/d" ); fMatchTree->Branch("deltaY_Trigger", &fDeltaY_Trigger, "deltaY_Trigger/d" ); - fMatchTree->Branch("radius_Trigger", &fRadius_Trigger, "dadius_Trigger/d" ); + fMatchTree->Branch("radius_Trigger", &fRadius_Trigger, "radius_Trigger/d" ); } //End MatchTree @@ -450,6 +457,7 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) // Cluster to Hit assns art::FindManyP cluster_hit_assns (clusterHandle, e, fPandoraLabel); + //For slice... for ( unsigned j = 0; j < nSlices; j++ ) { @@ -585,6 +593,11 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) //For flash... for ( int m = 0; m < nFlashes; m++ ) { auto & flash = flashVect[m]; + + double flashTime = flash->Time(); + + if(flashTimefFlashVetoWindow[1]) continue; + //Find index of flash that minimizes barycenter distance in XYZ place thisFlashCenterX = flash->XCenter(); thisFlashCenterY = flash->YCenter(); @@ -645,6 +658,7 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) if(chi2 < minChi2) { minChi2 = chi2; + fAngle = angle; matchIndex = idx; } } @@ -654,6 +668,8 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) const art::Ptr flashPtr { flashHandle, unsignedMatchIndex }; //Update match info updateFlashVars(flashPtr); + fChi2 = minChi2; + fScore = 1./fChi2; updateMatchInfo(sliceMatchInfo); sliceMatchInfoVector.push_back(sliceMatchInfo); art::Ptr const infoPtr = makeInfoPtr(matchInfoVector->size()); @@ -710,6 +726,9 @@ void TPCPMTBarycenterMatchProducer::InitializeSlice() { fDeltaY = -9999.; fDeltaZ = -9999.; fRadius = -9999.; + fChi2 = -9999.; + fScore = -99999.; + fAngle = -9999.; fDeltaZ_Trigger = -9999.; fDeltaY_Trigger = -9999.; fRadius_Trigger = -9999.; @@ -767,6 +786,9 @@ void TPCPMTBarycenterMatchProducer::updateMatchInfo(sbn::TPCPMTBarycenterMatch& matchInfo.deltaZ = fDeltaZ; matchInfo.radius = fRadius; matchInfo.radius_Trigger = fRadius_Trigger; + matchInfo.chi2 = fChi2; + matchInfo.score = fScore; + matchInfo.angle = fAngle; } //End updateMatchInfo() void TPCPMTBarycenterMatchProducer::GetPCA(std::vector const&x, std::vector const&y, std::vector const&weight, std::vector & PCA ) { From 67ee78c508c5fc32cfea39f0e71cb5dc97305f2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20S=C3=A1nchez=20Castillo?= <121103809+asanchezcastillo@users.noreply.github.com> Date: Fri, 7 Nov 2025 16:33:24 +0100 Subject: [PATCH 8/8] run bfm in one-to-many mode --- .../job/sbnd_tpcpmt3dbarycentermatching_config.fcl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl b/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl index d6f622c8a..4291fcb41 100644 --- a/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl +++ b/sbndcode/TPCPMTBarycenterMatching/job/sbnd_tpcpmt3dbarycentermatching_config.fcl @@ -23,7 +23,7 @@ TPCPMTBarycenterMatchProducer: YError: 23 // cm ZError: 23 // cm AngleError: 27 // deg - FlashVetoWindow: [-1500000,1500000] // ns + FlashVetoWindow: [-5,5] // ns module_type: "TPCPMTBarycenterMatchProducer" } @@ -31,4 +31,4 @@ TPCPMTBarycenterMatchProducerSCE: @local::TPCPMTBarycenterMatchProducer TPCPMTBarycenterMatchProducerSCE.PandoraLabel: "pandoraSCE" -END_PROLOG \ No newline at end of file +END_PROLOG