diff --git a/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt b/sbndcode/TPCPMTBarycenterMatching/CMakeLists.txt index 7158397dd..7ade54740 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,9 @@ set( MODULE_LIBRARIES larcoreobj::SummaryData larana::OpticalDetector_OpHitFinder larreco::Calorimetry + larsim::PhotonPropagation + larsim::PhotonPropagation_PhotonVisibilityService_service + larsim::OpticalPath lardata::Utilities nusimdata::SimulationBase nurandom::RandomUtils_NuRandomService_service @@ -31,6 +35,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..70fea3aa7 100644 --- a/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc +++ b/sbndcode/TPCPMTBarycenterMatching/TPCPMTBarycenterMatching_module.cc @@ -20,6 +20,7 @@ #include "art/Framework/Principal/Handle.h" #include "art/Framework/Principal/Run.h" #include "art/Framework/Principal/SubRun.h" +#include "art/Utilities/make_tool.h" #include "canvas/Utilities/InputTag.h" #include "fhiclcpp/ParameterSet.h" #include "messagefacility/MessageLogger/MessageLogger.h" @@ -39,7 +40,11 @@ #include "larcore/Geometry/Geometry.h" #include "larcore/Geometry/WireReadout.h" #include "larcorealg/Geometry/GeometryCore.h" +#include "larsim/PhotonPropagation/SemiAnalyticalModel.h" +#include "larsim/PhotonPropagation/OpticalPathTools/OpticalPath.h" #include "larcore/CoreUtils/ServiceUtil.h" +#include "sbndcode/OpDetSim/sbndPDMapAlg.hh" + //Data product includes #include "lardataobj/RecoBase/OpHit.h" @@ -48,13 +53,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 +218,31 @@ 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 @@ -232,6 +258,7 @@ 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) @@ -239,16 +266,32 @@ 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 fDeltaX_Trigger; ///< | Triggering flash X center - charge X center | (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) - 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; + std::shared_ptr _optical_path_tool; }; @@ -259,9 +302,19 @@ TPCPMTBarycenterMatchProducer::TPCPMTBarycenterMatchProducer(fhicl::ParameterSet fOpFlashesModuleLabel(p.get>("OpFlashesModuleLabel")), fPandoraLabel(p.get("PandoraLabel")), fCollectionOnly(p.get("CollectionOnly", true)), - fVerbose(p.get("Verbose", false)), - fFillMatchTree(p.get("FillMatchTree", false)), - fDo3DMatching(p.get("Do3DMatching",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 { // Call appropriate produces<>() functions here. @@ -274,14 +327,15 @@ 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 @@ -305,17 +359,36 @@ 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("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 + //Fill the OpDet positions + for(unsigned int opch=0; opch("VUVHits"); + _vis_params = p.get("VIVHits"); + _optical_path_tool = std::shared_ptr(art::make_tool(p.get("OpticalPathTool"))); + _semi_model = std::make_unique(_vuv_params, _vis_params, _optical_path_tool, true, false); } void TPCPMTBarycenterMatchProducer::produce(art::Event& e) @@ -324,93 +397,148 @@ 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); - art::FindManyP fmTPCHits(sliceHandle, e, fPandoraLabel); - art::FindManyP fmPFPs(sliceHandle, e, fPandoraLabel); + //Fetch slices, TPC hits, and PFPs; pointer vector needed for generating associations + ::art::Handle> sliceHandle; + e.getByLabel(fPandoraLabel, sliceHandle); + + // Slice to PFP assns 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 }; + //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 ( 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); - int nHits = tpcHitsVec.size(); - int nPFPs = pfpsVec.size(); + std::vector hit_z; + std::vector hit_y; + std::vector hit_weight; - //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.; 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 @@ -419,6 +547,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 +563,61 @@ 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; - 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(); + //For flash... for ( int m = 0; m < nFlashes; m++ ) { - const recob::OpFlash &flash = (*flashHandle).at(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(); - 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()); @@ -476,7 +683,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,6 +693,7 @@ void TPCPMTBarycenterMatchProducer::produce(art::Event& e) minNPes = nPEs; } } + if(maxFlashIdx!=-1) { sliceAssns->addSingle(infoPtrVector[maxFlashIdx], slicePtr); @@ -518,11 +726,12 @@ void TPCPMTBarycenterMatchProducer::InitializeSlice() { fFlashWidthY = -9999.; fFlashWidthZ = -9999.; fDeltaT = -9999.; - fDeltaX = -9999.; fDeltaY = -9999.; fDeltaZ = -9999.; fRadius = -9999.; - fDeltaX_Trigger = -9999.; + fChi2 = -9999.; + fScore = -99999.; + fAngle = -9999.; fDeltaZ_Trigger = -9999.; fDeltaY_Trigger = -9999.; fRadius_Trigger = -9999.; @@ -538,14 +747,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 +772,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() @@ -592,7 +789,119 @@ 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