From 47efd4878c86dce841b5b6b38b2a77f68d2f0043 Mon Sep 17 00:00:00 2001 From: nathanielerowe <70993723+nathanielerowe@users.noreply.github.com> Date: Tue, 16 Sep 2025 07:53:26 -0500 Subject: [PATCH] Revert "Feature/acastill tpcpmt bfm upgrade" --- .../TPCPMTBarycenterMatching/CMakeLists.txt | 4 - .../TPCPMTBarycenterMatching_module.cc | 484 ++++-------------- ...sbnd_tpcpmt3dbarycentermatching_config.fcl | 31 +- 3 files changed, 100 insertions(+), 419 deletions(-) diff --git a/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt b/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt index 9a409d8f0..7158397dd 100644 --- a/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt +++ b/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt @@ -1,6 +1,5 @@ set( MODULE_LIBRARIES sbnobj::Common_Trigger - sbndcode_OpDetSim larcorealg::Geometry larcore::Geometry_Geometry_service lardataobj::RecoBase @@ -8,8 +7,6 @@ set( MODULE_LIBRARIES larcoreobj::SummaryData larana::OpticalDetector_OpHitFinder larreco::Calorimetry - larsim::PhotonPropagation - larsim::PhotonPropagation_PhotonVisibilityService_service lardata::Utilities nusimdata::SimulationBase nurandom::RandomUtils_NuRandomService_service @@ -34,7 +31,6 @@ 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 a9f260fe3..8da5761d5 100644 --- a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc +++ b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc @@ -39,10 +39,7 @@ #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" @@ -51,19 +48,13 @@ #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 @@ -216,31 +207,16 @@ 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 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 - 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 - 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; + // Event-level data members int fRun; ///< Number of the run being processed int fEvent; ///< Number of the event being processed @@ -256,7 +232,6 @@ 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) @@ -264,31 +239,16 @@ class TPCPMTBarycenterMatchProducer : public art::EDProducer { 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 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 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 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) + TTree* fMatchTree; ///< Tree to store all match 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; }; @@ -299,19 +259,9 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet fOpFlashesModuleLabel(p.get>("OpFlashesModuleLabel")), fPandoraLabel(p.get("PandoraLabel")), 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")), - 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 + fVerbose(p.get("Verbose", false)), + fFillMatchTree(p.get("FillMatchTree", false)), + fDo3DMatching(p.get("Do3DMatching",true)) { // Call appropriate produces<>() functions here. @@ -324,15 +274,14 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet << "Module has been missconfigured. Number of OpFlashesModuleLabel must be 2!"; } - art::ServiceHandle tfs; if ( fFillMatchTree ) { - art::ServiceHandle tfs; - fMatchTree = tfs->make("matchTree","TPC Slice - OpFlash Matching Analysis"); + art::ServiceHandle tfs; + fMatchTree = tfs->make("matchTree","TPC Slice - OpFlash Matching Analysis"); //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 @@ -356,36 +305,17 @@ 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("angle", &fAngle, "angle/d" ); - fMatchTree->Branch("chi2", &fChi2, "chi2/d" ); - fMatchTree->Branch("score", &fScore, "score/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, "radius_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) @@ -394,148 +324,93 @@ 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 - */ - - //Fetch slices, TPC hits, and PFPs; pointer vector needed for generating associations - ::art::Handle> sliceHandle; - e.getByLabel(fPandoraLabel, sliceHandle); +/* ~~~~~~~~~~~~~~~~~~~~ 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 + */ - // Slice to PFP assns + //Fetch slices, TPC hits, and PFPs; pointer vector needed for generating associations + art::Handle const sliceHandle + = e.getHandle>(fPandoraLabel); + art::FindManyP fmTPCHits(sliceHandle, e, fPandoraLabel); + art::FindManyP fmPFPs(sliceHandle, e, fPandoraLabel); unsigned nSlices = (*sliceHandle).size(); - ::art::Handle> pfpHandle; - e.getByLabel(fPandoraLabel, pfpHandle); - - // PFP to metadata assns - art::FindManyP slice_pfp_assns(sliceHandle, e, fPandoraLabel); - - // Hit to PFP assns - 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); - - //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 }; - //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 = 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 each TPC + std::vector sliceMatchInfoVector; + std::vector> infoPtrVector; + std::vector> flashPtrVector; + fSliceNum = j; + const art::Ptr slicePtr { sliceHandle, j }; for ( size_t tpc=0; tpc<2 ; tpc ++) { fTPC = tpc; InitializeSlice(); sbn::TPCPMTBarycenterMatch sliceMatchInfo; updateMatchInfo(sliceMatchInfo); - art::FindOne f1SpacePoint(tpcHitsVec, e, fPandoraLabel); - std::vector hit_z; - std::vector hit_y; - std::vector hit_weight; + const std::vector> &tpcHitsVec = fmTPCHits.at(j); + const std::vector> &pfpsVec = fmPFPs.at(j); + art::FindOne f1SpacePoint(tpcHitsVec, e, fPandoraLabel); 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; + } + } double thisCharge; double sumCharge = 0.; TVector3 sumPos {0.,0.,0.}; TVector3 sumPosSqr {0.,0.,0.}; - size_t maxChargePlaneIdx=99999; - - // 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++ ) { - 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 @@ -544,10 +419,6 @@ 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... @@ -560,61 +431,27 @@ 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); - //Get the visibility map of the charge barycenter - geo::Point_t SliceXYZ = {fChargeCenterX, fChargeCenterY, fChargeCenterZ}; - 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; - // --- 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(); - + art::Handle flashHandle = e.getHandle>(fOpFlashesModuleLabel[tpc]); + int nFlashes = (*flashHandle).size(); //For flash... for ( int m = 0; m < nFlashes; m++ ) { - auto & flash = flashVect[m]; - - double flashTime = flash->Time(); - - if(flashTimefFlashVetoWindow[1]) continue; - + const recob::OpFlash &flash = (*flashHandle).at(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) ); - if ( thisDistance < minDistance ) { minDistance = thisDistance; matchIndex = m; } - // Fill the vector with the idxs of the flashes - if(thisDistance 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); - - double FlashLight = this->GetFlashLight(flash->TotalPE(), direct_visibility, reflect_visibility); - fFlashTimeDebug = flash->Time(); - double lightChargeRatio = (FlashLight/sliceCharge); - std::vector LightPCA = {0.,0.}; - GetPCA(ophit_z, ophit_y, ophit_weight, LightPCA); - 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 - 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(lightChargeRatio < fLightChargeRatioBounds[0] || lightChargeRatio > fLightChargeRatioBounds[1]) continue; - if(chi2 < minChi2) - { - minChi2 = chi2; - fAngle = angle; - matchIndex = idx; - } - } - //Best match flash pointer unsigned unsignedMatchIndex = matchIndex; 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()); @@ -680,7 +476,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(); @@ -690,7 +486,6 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) minNPes = nPEs; } } - if(maxFlashIdx!=-1) { sliceAssns->addSingle(infoPtrVector[maxFlashIdx], slicePtr); @@ -723,12 +518,11 @@ void TPCPMTBarycenterMatchProducer::InitializeSlice() { fFlashWidthY = -9999.; fFlashWidthZ = -9999.; fDeltaT = -9999.; + fDeltaX = -9999.; fDeltaY = -9999.; fDeltaZ = -9999.; fRadius = -9999.; - fChi2 = -9999.; - fScore = -99999.; - fAngle = -9999.; + fDeltaX_Trigger = -9999.; fDeltaZ_Trigger = -9999.; fDeltaY_Trigger = -9999.; fRadius_Trigger = -9999.; @@ -744,8 +538,14 @@ void TPCPMTBarycenterMatchProducer::updateChargeVars(double sumCharge, TVector3 fChargeWidthZ = std::sqrt( sumPosSqr[2]/sumCharge - (sumPos[2]/sumCharge)*(sumPos[2]/sumCharge) ); fChargeTotal = sumCharge; if ( triggerFlashCenter[1] != -9999. ) { - fDeltaY_Trigger = abs(triggerFlashCenter[0] - fChargeCenterY); - fDeltaZ_Trigger = abs(triggerFlashCenter[1] - fChargeCenterZ); + 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 fRadius_Trigger = std::hypot(fDeltaY_Trigger, fDeltaZ_Trigger); } } //End updateChargeVars() @@ -769,6 +569,12 @@ 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() @@ -786,119 +592,7 @@ 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 ) { - - 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]; - } - - //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) { - Eigen::Vector2d pt = centered_points.row(i); - cov += std::abs(weight[i]) * (pt * pt.transpose()); - } - cov /= weight_sum; - - // Get eigenvalues / eigenvectors - Eigen::SelfAdjointEigenSolver solver(cov); - Eigen::VectorXd eigvals = solver.eigenvalues(); - Eigen::Matrix2d eigvecs = solver.eigenvectors(); - - // Get Principal component (largest eigenvalue) - 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 bestHits = std::max_element(plane_hits.begin(), plane_hits.end()) - plane_hits.begin(); - double _comp_charge = plane_charge.at(bestHits); - - 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