From 7f7a897f0b7c26c22885b4c62e15c3d8ea99c426 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 09:44:50 -0600 Subject: [PATCH 01/41] Adding blip data types to sbnobj to help with CAF addition --- sbnobj/SBND/Blip/BlipDataTypes.h | 171 +++++++++++++++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100644 sbnobj/SBND/Blip/BlipDataTypes.h diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h new file mode 100644 index 00000000..79d7f8dc --- /dev/null +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -0,0 +1,171 @@ +#include "lardataobj/RecoBase/Hit.h" +#include "nusimdata/SimulationBase/MCParticle.h" + +#include + +typedef std::vector vint_t; +typedef std::vector vbool_t; +typedef std::vector vfloat_t; +typedef std::set si_t; +typedef std::map mif_t; + +const int kNplanes = 3; +const int kNTPCs = 2; + +namespace blip { + + //################################################### + // Data structures + //################################################### + + struct ParticleInfo { + simb::MCParticle particle; + int trackId = -9; + int index = -9; + int isPrimary = -9; + int numTrajPts = -9; + double depEnergy = -9; + int depElectrons = -9; + double numElectrons = -9; + double mass = -9; + double E = -9; + double endE = -9; + double KE = -9; + double endKE = -9; + double P = -9; + double Px = -9; + double Py = -9; + double Pz = -9; + double pathLength = -9; + double time = -9; + double endtime = -9; + TVector3 startPoint; + TVector3 endPoint; + TVector3 position; + }; + + // True energy depositions + struct TrueBlip { + int ID = -9; // unique blip ID + int Cryostat = -9; // Cryostat ID + int TPC = -9; // TPC ID + float Time = -9; // time of particle interaction + int TimeTick = -9; // time tick + float DriftTime = -9; // drift time [us] + float Energy = 0; // energy dep [MeV] + int DepElectrons = 0; // deposited electrons + int NumElectrons = 0; // electrons reaching wires + int LeadG4ID = -9; // lead G4 track ID + int LeadG4Index = -9; // lead G4 track index + int LeadG4PDG = -9; // lead G4 PDG + float LeadCharge = -9; // lead G4 charge dep + TVector3 Position; // XYZ position + mif_t G4ChargeMap; + mif_t G4PDGMap; + }; + + struct HitInfo { + int hitid = -9; + int cryo = -9; + int tpc = -9; + int plane = -9; + int wire = -9; + int chan = -9; + float amp = -9; + float rms = -9; + int trkid = -9; + int shwrid = -9; + int clustid = -9; + int blipid = -9; + bool ismatch = false; + float integralADC = -999; // [ADCs] from integral + float sigmaintegral = -999; + float sumADC = -999; // [ADCs] from sum + float charge = -999; // [e-] + float peakTime = -999999; + float driftTime = -999999; // [tick] + float gof = -9; + int g4trkid = -9; + int g4pdg = -999; + int g4charge = -999; // [e-] + float g4frac = -99; + float g4energy = -999; // [MeV] + }; + + struct HitClust { + int ID = -9; + bool isValid = false; + int CenterChan = -999; + int CenterWire = -999; + bool isTruthMatched = false; + bool isMerged = false; + bool isMatched = false; + int DeadWireSep = 99; + int Cryostat = -9; + int TPC = -9; + int Plane = -9; + int NHits = -9; + int NWires = -9; + float ADCs = -999; + float Amplitude = -999; + float Charge = -999; + float SigmaCharge = -999; + float TimeTick = -999; + float Time = -999; + float StartHitTime = -999; + float EndHitTime = -999; + float StartTime = -999; + float EndTime = -999; + float Timespan = -999; + float RMS = -999; + int StartWire = -999; + int EndWire = -999; + int NPulseTrainHits = -9; + float GoodnessOfFit = -999; + int BlipID = -9; + int EdepID = -9; + si_t HitIDs; + si_t Wires; + si_t Chans; + si_t G4IDs; + + std::map IntersectLocations; + }; + + struct Blip { + + int ID = -9; // Blip ID / index + bool isValid = false; // Blip passes basic checks + int Cryostat = -9; // Cryostat + int TPC = -9; // TPC + int NPlanes = -9; // Num. matched planes + int MaxWireSpan = -9; // Maximum span of wires on any plane cluster + float TimeTick = -999; // Readout time [ticks] + float Time = -999; // Drift time [us] + float Charge = -9; // Charge on calorimetry plane + float Energy = -999; // Energy (const dE/dx, fcl-configurable) + float EnergyESTAR = -999; // Energy (ESTAR method from ArgoNeuT) + float EnergyPSTAR = -999; // Energy (PSTAR method similar with ESTAR method from ArgoNeuT) + float ProxTrkDist = -9; // Distance to cloest track + int ProxTrkID = -9; // ID of closest track + bool inCylinder = false; // Is it in a cone/cylinder region? + + TVector3 Position; // 3D position TVector3 + float SigmaYZ = -9.; // Uncertainty in YZ intersect [cm] + float dX = -9; // Equivalent length along drift direction [cm] + float dYZ = -9; // Approximate length scale in YZ space [cm] + + // Plane/cluster-specific information + blip::HitClust clusters[kNplanes]; + + // Truth-matched energy deposition + blip::TrueBlip truth; + + // Prototype getter functions + double X() { return Position.X(); } + double Y() { return Position.Y(); } + double Z() { return Position.Z(); } + + }; + +} From 18d8269518ddeddeb9ad902d2daa091b05cd276e Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 09:55:48 -0600 Subject: [PATCH 02/41] Moving blip utils over --- sbnobj/SBND/Blip/BlipUtils.cc | 749 ++++++++++++++++++ sbnobj/SBND/Blip/BlipUtils.h | 92 +++ sbnobj/SBND/Blip/CMakeLists.txt | 45 ++ .../Blip/{BlipDataTypes.h => DataTypes.h} | 3 + sbnobj/SBND/Blip/classes.h | 36 + sbnobj/SBND/Blip/classes_def.xml | 14 + sbnobj/SBND/CMakeLists.txt | 1 + 7 files changed, 940 insertions(+) create mode 100644 sbnobj/SBND/Blip/BlipUtils.cc create mode 100644 sbnobj/SBND/Blip/BlipUtils.h create mode 100644 sbnobj/SBND/Blip/CMakeLists.txt rename sbnobj/SBND/Blip/{BlipDataTypes.h => DataTypes.h} (99%) create mode 100644 sbnobj/SBND/Blip/classes.h create mode 100644 sbnobj/SBND/Blip/classes_def.xml diff --git a/sbnobj/SBND/Blip/BlipUtils.cc b/sbnobj/SBND/Blip/BlipUtils.cc new file mode 100644 index 00000000..1b1e674f --- /dev/null +++ b/sbnobj/SBND/Blip/BlipUtils.cc @@ -0,0 +1,749 @@ +#include "BlipUtils.h" + +namespace BlipUtils { + + //============================================================================ + // Find total visible energy deposited in the LAr, and number of electrons deposited + // and drifted to the anode. + /* + void CalcTotalDep(float& energy, int& ne_dep, float& ne_anode, SEDVec_t& sedvec){ + + // energy and electrons deposited + energy = 0; + ne_dep = 0; + for(auto& sed : sedvec ) { + energy += sed->Energy(); + ne_dep += sed->NumElectrons(); + } + + // electrons drifted to collection plane wires + art::ServiceHandle geom; + ne_anode = 0; + for(auto const &chan : art::ServiceHandle()->SimChannels()) { + if( geom->View(chan->Channel()) != geo::kW ) continue; + for(auto const &tdcide : chan->TDCIDEMap() ) { + for(const auto& ide : tdcide.second) ne_anode += ide.numElectrons; + } + } + } + */ + + + //=========================================================================== + // Provided a MCParticle, calculate everything we'll need for later calculations + // and save into ParticleInfo object + void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo, SEDVec_t& sedvec, int caloPlane){ + + // Get important info and do conversions + pinfo.particle = part; + pinfo.trackId = part.TrackId(); + pinfo.isPrimary = (int)(part.Process() == "primary"); + pinfo.mass = /*GeV->MeV*/1e3 * part.Mass(); + pinfo.E = /*GeV->MeV*/1e3 * part.E(); + pinfo.endE = /*GeV->MeV*/1e3 * part.EndE(); + pinfo.KE = /*GeV->MeV*/1e3 * (part.E()-part.Mass()); + pinfo.endKE = /*GeV->MeV*/1e3 * (part.EndE()-part.Mass()); + pinfo.P = /*GeV->MeV*/1e3 * part.Momentum().Vect().Mag(); + pinfo.Px = /*GeV->MeV*/1e3 * part.Px(); + pinfo.Py = /*GeV->MeV*/1e3 * part.Py(); + pinfo.Pz = /*GeV->MeV*/1e3 * part.Pz(); + pinfo.time = /*ns ->mus*/1e-3 * part.T(); + pinfo.endtime = /*ns ->mus*/1e-3 * part.EndT(); + pinfo.numTrajPts = part.NumberTrajectoryPoints(); + + // Pathlength (in AV) and start/end point + pinfo.pathLength = PathLength( part, pinfo.startPoint, pinfo.endPoint); + + // Central position of trajectory + pinfo.position = 0.5*(pinfo.startPoint+pinfo.endPoint); + + // Energy/charge deposited by this particle, found using SimEnergyDeposits + pinfo.depEnergy = 0; + pinfo.depElectrons = 0; + for(auto& sed : sedvec ) { + if( sed->TrackID() == part.TrackId() ) { + pinfo.depEnergy += sed->Energy(); + pinfo.depElectrons += sed->NumElectrons(); + } + } + + return; + + } + + //=================================================================== + // Provided a vector of all particle information for event, fill a + // vector of true blips + void MakeTrueBlips( std::vector& pinfo, std::vector& trueblips ) { + + art::ServiceHandle geom; + auto const& detProp = art::ServiceHandle()->DataForJob(); + auto const& clockData = art::ServiceHandle()->DataForJob(); + + for(size_t i=0; i0) ? tb.Time/clockData.TPCClock().TickPeriod() : 0; + auto point = geo::Point_t{tb.Position.X(),tb.Position.Y(),tb.Position.Z()}; + auto const& tpcID = geom->FindTPCAtPosition(point); + auto const& planeID = geo::PlaneID{tpcID, 0}; + float tick_calc = (float)detProp.ConvertXToTicks(tb.Position.X(),planeID); + tb.DriftTime = tick_calc*clockData.TPCClock().TickPeriod() + clockData.TriggerOffsetTPC(); + + tb.ID = trueblips.size(); + trueblips.push_back(tb); + + } + + } + + + //==================================================================== + void GrowTrueBlip( blip::ParticleInfo& pinfo, blip::TrueBlip& tblip ) { + + simb::MCParticle& part = pinfo.particle; + + // Skip neutrons, photons + if( part.PdgCode() == 2112 || part.PdgCode() == 22 ) return; + + // Check that path length isn't zero + if( !pinfo.pathLength ) return; + + // If this is a new blip, initialize + if( !tblip.G4ChargeMap.size() ) { + tblip.Position = pinfo.position; + tblip.Time = pinfo.time; + + // .. otherwise, check that the new particle + // creation time is comparable to existing blip. + // then calculate new energy-weighted position. + } else if ( fabs(tblip.Time-pinfo.time) < 3 ) { + float totE = tblip.Energy + pinfo.depEnergy; + float w1 = tblip.Energy/totE; + float w2 = pinfo.depEnergy/totE; + tblip.Position = w1*tblip.Position + w2*pinfo.position; + tblip.Time = w1*tblip.Time + w2*pinfo.time; + tblip.LeadCharge = pinfo.depElectrons; + // ... if the particle isn't a match, show's over + } else { + return; + } + + tblip.Energy += pinfo.depEnergy; + tblip.DepElectrons+= pinfo.depElectrons; + tblip.NumElectrons+= std::max(0.,pinfo.numElectrons); + + tblip.G4ChargeMap[part.TrackId()] += pinfo.depElectrons; + tblip.G4PDGMap[part.TrackId()] = part.PdgCode(); + if(pinfo.depElectrons > tblip.LeadCharge ) { + tblip.LeadCharge = pinfo.depElectrons; + tblip.LeadG4Index = pinfo.index; + tblip.LeadG4ID = part.TrackId(); + tblip.LeadG4PDG = part.PdgCode(); + } + } + + + //==================================================================== + // Merge blips that are close + void MergeTrueBlips(std::vector& vtb, float dmin){ + if( dmin <= 0 ) return; + std::vector vtb_merged; + std::vector isGrouped(vtb.size(),false); + + for(size_t i=0; i 5 ) continue; + float d = (blip_i.Position-blip_j.Position).Mag(); + if( d < dmin ) { + isGrouped.at(j) = true; + //float totE = blip_i.Energy + blip_j.Energy; + float totQ = blip_i.DepElectrons + blip_j.DepElectrons; + float w1 = blip_i.DepElectrons/totQ; + float w2 = blip_j.DepElectrons/totQ; + blip_i.Energy += blip_j.Energy; + blip_i.Position = w1*blip_i.Position + w2*blip_j.Position; + blip_i.DriftTime = w1*blip_i.DriftTime+ w2*blip_j.DriftTime; + blip_i.Time = w1*blip_i.Time + w2*blip_j.Time; + blip_i.DepElectrons += blip_j.DepElectrons; + if( blip_j.NumElectrons ) blip_i.NumElectrons += blip_j.NumElectrons; + + blip_i.G4ChargeMap.insert(blip_j.G4ChargeMap.begin(), blip_j.G4ChargeMap.end()); + blip_i.G4PDGMap.insert(blip_j.G4PDGMap.begin(), blip_j.G4PDGMap.end()); + + if( blip_j.LeadCharge > blip_i.LeadCharge ) { + blip_i.LeadCharge = blip_j.LeadCharge; + blip_i.LeadG4ID = blip_j.LeadG4ID; + blip_i.LeadG4Index = blip_j.LeadG4Index; + blip_i.LeadG4PDG = blip_j.LeadG4PDG; + } + }//d < dmin + }//loop over blip_j + blip_i.ID = vtb_merged.size(); + vtb_merged.push_back(blip_i); + } + vtb.clear(); + vtb = vtb_merged; + } + + + //================================================================= + blip::HitClust MakeHitClust(std::vector const& hitinfoVec){ + + blip::HitClust hc; + if( hitinfoVec.size() ) { + int tpc = hitinfoVec[0].tpc; + int plane = hitinfoVec[0].plane; + int cryo = hitinfoVec[0].cryo; + + // check that all hits are on same plane; + for(auto& h : hitinfoVec ) { + if( h.cryo != cryo ) return hc; + if( h.tpc != tpc ) return hc; + if( h.plane != plane ) return hc; + } + + // initialize values + hc.Cryostat = cryo; + hc.TPC = tpc; + hc.Plane = plane; + hc.ADCs = 0; + hc.Charge = 0; + hc.SigmaCharge = 0; + hc.Amplitude = 0; + hc.NPulseTrainHits = 0; + float startTime = 9e9; + float endTime = -9e9; + float weightedTick = 0; + float weightedTime = 0; + float weightedGOF = 0; + //float weightedRatio = 0; + float qGOF = 0; + + // store hit times, charges, and RMS + std::vector tvec; + std::vector qvec; + std::vector dqvec; + std::vector rmsvec; + + // grow our hit cluster! + for(auto& hitinfo : hitinfoVec ) { + if( hc.HitIDs.find(hitinfo.hitid) != hc.HitIDs.end() ) continue; + hc.HitIDs .insert(hitinfo.hitid); + hc.Wires .insert(hitinfo.wire); + hc.Chans .insert(hitinfo.chan); + float q = (hitinfo.charge > 0)? hitinfo.charge : 0; + float integral = hitinfo.integralADC; + float sigma = hitinfo.sigmaintegral; + float dq = (integral != 0 && sigma>0)? (sigma/integral)*q : 0; + hc.Charge += q; + hc.ADCs += hitinfo.integralADC; + hc.Amplitude = std::max(hc.Amplitude, hitinfo.amp ); + weightedTick += q*hitinfo.peakTime; + weightedTime += q*hitinfo.driftTime; + startTime = std::min(startTime, hitinfo.driftTime-hitinfo.rms); + endTime = std::max(endTime, hitinfo.driftTime+hitinfo.rms); + tvec .push_back(hitinfo.driftTime); + qvec .push_back(q); + dqvec .push_back(dq); + rmsvec .push_back(hitinfo.rms); + if( hitinfo.g4trkid >= 0 ) hc.G4IDs.insert(hitinfo.g4trkid); + if( hitinfo.gof < 0 ) { + hc.NPulseTrainHits++; + } else { + weightedGOF += q*hitinfo.gof; + qGOF += q; + } + }//endloop over hits + + // mean goodness of fit + if( qGOF ) hc.GoodnessOfFit = weightedGOF/qGOF; + + // calculate other quantities + hc.NHits = hc.HitIDs.size(); + hc.NWires = hc.Wires.size(); + hc.CenterWire =(*hc.Wires.begin()+*hc.Wires.rbegin())/2.; + hc.CenterChan =(*hc.Chans.begin()+*hc.Chans.rbegin())/2.; + hc.StartWire = *hc.Wires.begin(); + hc.EndWire = *hc.Wires.rbegin(); + hc.StartTime = startTime; + hc.EndTime = endTime; + hc.Timespan = hc.EndTime - hc.StartTime; + hc.Time = weightedTime / hc.Charge; + hc.TimeTick = weightedTick / hc.Charge; + + // overall cluster RMS and uncertainty in charge + float sig_sumSq = 0; + float dt_sumSq = 0; + float dq = 0; + for(size_t i=0; i 0 hits + + // mark the cluster as valid and ship it out + hc.isValid = true; + return hc; + } + + + //================================================================= + blip::Blip MakeBlip( std::vector const& hcs, + detinfo::DetectorPropertiesData const& detProp, + detinfo::DetectorClocksData const& clockData ){ + + art::ServiceHandle wireReadoutGeom; + blip::Blip newblip; + + // ------------------------------------------------ + // Must be 1-3 clusts (no more, no less!) + if( hcs.size() > 3 || hcs.size() < 1 ) return newblip; + + // ------------------------------------------------ + // All hits must be in same TPC, and no 2 planes can be repeated + std::set planeIDs; + for(size_t i=0; i()->DataForJob(); + //auto const& clockData = art::ServiceHandle()->DataForJob(); + float driftVelocity = detProp.DriftVelocity(detProp.Efield(0),detProp.Temperature()); + float tick_to_cm = clockData.TPCClock().TickPeriod() * driftVelocity; + + + newblip.Cryostat = hcs[0].Cryostat; + newblip.TPC = hcs[0].TPC; + newblip.NPlanes = planeIDs.size(); + + // ------------------------------------------------ + /// Look for valid wire intersections between + // central-most hits in each cluster + std::vector wirex; + for(size_t i=0; iGet().Plane(geo::PlaneID{(unsigned int)hcs[i].Cryostat, (unsigned int)hcs[i].TPC, (unsigned int)hcs[i].Plane}); + double wirepitch = planegeo.WirePitch(); + // use view with the maximal wire extent to calculate transverse (YZ) length + if( hcs[i].NWires > newblip.MaxWireSpan ) { + newblip.MaxWireSpan = hcs[i].NWires; + newblip.dYZ = hcs[i].NWires * wirepitch; + } + + for(size_t j=i+1; jsecond.Y()); + intsec_p.SetZ(hcs[i].IntersectLocations.find(hcs[j].ID)->second.Z()); + } else { + std::vector i_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcs[i].CenterChan); + std::vector j_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcs[j].CenterChan); + + match3d = wireReadoutGeom->Get().WireIDsIntersect(i_wireids.at(0), j_wireids.at(0), intsec_p); + } + + if( match3d ) { + TVector3 a(0., intsec_p.Y(), intsec_p.Z()); + wirex.push_back(a); + newblip.clusters[pli] = hcs[i]; + newblip.clusters[plj] = hcs[j]; + } + } + } + + // Require some number of intersection points. + if( !wirex.size() ) return newblip; + + // Loop over the intersection points and calculate average position in + // YZ-plane, as well as the mean difference between intersection points. + newblip.Position.SetXYZ(0,0,0); + if( wirex.size() == 1 ) { + newblip.Position= wirex[0]; + } else { + newblip.SigmaYZ = 0; + double fact = 1./wirex.size(); + for(auto& v : wirex ) newblip.Position += v * fact; + for(auto& v : wirex ) newblip.SigmaYZ += (v-newblip.Position).Mag() * fact; + // Ensure that difference between intersection points is + // consistent with the maximal wire extent + if( newblip.SigmaYZ > std::max(1.,0.5*newblip.dYZ) ) return newblip; + } + + // Calculate mean drift time and X-position + // (note that the 'time' of each of the hit clusters + // have already been corrected for plane-to-plane offsets) + newblip.TimeTick= 0; + newblip.Time = 0; + newblip.dX = 0; + float vsize = (float)hcs.size(); + for(auto hc : hcs ) { + newblip.TimeTick += hc.TimeTick / vsize; + newblip.Time += hc.Time / vsize; + newblip.dX = std::max((float)(hc.EndTime-hc.StartTime)*tick_to_cm, newblip.dX); + } + + //auto const& tpcID = geo::TPCID(geo::CryostatID(newblip.Cryostat),newblip.TPC); + //auto const& planeID = art::ServiceHandle()->GetBeginPlaneID(tpcID); + //newblip.Position .SetX(detProp.ConvertTicksToX(newblip.TimeTick, planeID)); + + // convert ticks to X + auto const& cryostat= art::ServiceHandle()->Cryostat(geo::CryostatID(newblip.Cryostat)); + auto const& tpcgeom = cryostat.TPC(newblip.TPC); + auto tpcID = tpcgeom.ID(); + auto const& planegeo = wireReadoutGeom->Get().Plane(geo::PlaneID{tpcID, 0}); + auto const xyz = planegeo.GetCenter(); + int dirx = DriftDirX(tpcgeom); + + newblip.Position.SetX( xyz.X() + dirx * tick_to_cm * newblip.Time ); + + // this should ALREADY be accounted for at the hit-processing level in BlipRecoAlg, + // through the use of GetXTicksOffset... + //float offset_ticks = clockData.TriggerOffsetTPC() / clockData.TPCClock().TickPeriod(); + //float driftTicks = newblip.Time + clockData.TriggerOffsetTPC(); + //std::cout<<"blip time "< pi_serv; + const sim::ParticleList& plist = pi_serv->ParticleList(); + if( particleID == ancestorID ) return true; + if( particleID < ancestorID ) return false; + if( !plist.HasParticle(ancestorID) ) return false; + while( particleID > ancestorID ) { + + simb::MCParticle p = pi_serv->TrackIdToParticle(particleID); + if( !plist.HasParticle(p.Mother() ) ) { return false; } + + simb::MCParticle pM = pi_serv->TrackIdToParticle(p.Mother()); + if ( pM.TrackId() == ancestorID ) { return true; } + else if ( breakAtPhots == true && pM.PdgCode() == 22 ) { return false; } + else if ( pM.Process() == "primary" || pM.TrackId() == 1 ) { return false; } + else if ( pM.Mother() == 0 ) { return false; } + else { particleID = pM.TrackId(); } + } + + return false; + } + + //=================================================================== + int DriftDirX(geo::TPCGeo const& tpcgeom) { + return ((tpcgeom.DriftSign() == geo::DriftSign::Negative) ? +1.0 : -1.0); + } + + //==================================================================== + bool DoHitsOverlap(art::Ptr const& hit1, art::Ptr const& hit2){ + if( hit1->WireID() != hit2->WireID() ) return false; + float t1 = hit1->PeakTime(); + float t2 = hit2->PeakTime(); + float sig = std::max(hit1->RMS(),hit2->RMS()); + if( fabs(t1-t2) < 1.0*sig ) return true; + else return false; + } + + + //==================================================================== + bool DoHitClustsOverlap(blip::HitClust const& hc1, blip::HitClust const& hc2){ + + // only match across different wires in same TPC + if( hc1.TPC != hc2.TPC ) return false; + + if( hc1.StartTime <= hc2.EndTime + && hc2.StartTime <= hc1.EndTime ) return true; + else return false; + } + bool DoHitClustsOverlap(blip::HitClust const& hc1, float t1, float t2 ){ + blip::HitClust hc2; + hc2.TPC = hc1.TPC; + hc2.StartTime = t1; + hc2.EndTime = t2; + return DoHitClustsOverlap(hc1,hc2); + } + + //==================================================================== + // Calculates the level of time overlap between two clusters + float CalcHitClustsOverlap(blip::HitClust const& hc1, blip::HitClust const& hc2){ + return CalcOverlap(hc1.StartTime,hc1.EndTime,hc2.StartTime,hc2.EndTime); + } + + float CalcOverlap(const float& x1, const float& x2, const float& y1, const float& y2){ + float full_range = std::max(x2,y2) - std::min(x1,y1); + float sum = (x2-x1) + (y2-y1); + float overlap = std::max(float(0), sum-full_range); + if( overlap > 0 ) return 2. * overlap / sum; + else return -1; + } + + //==================================================================== + bool DoHitClustsMatch(blip::HitClust const& hc1, blip::HitClust const& hc2, float minDiffTicks = 2){ + if( fabs(hc1.Time-hc2.Time) < minDiffTicks ) return true; + else return false; + } + + //==================================================================== + // This function calculates the leading MCParticle ID contributing to a hit and the + // fraction of that hit's energy coming from that particle. + /* + void HitTruth(art::Ptr const& hit, int& truthid, float& truthidEnergyFrac, float& energy,float& numElectrons){ + // Get associated sim::TrackIDEs for this hit + std::vector trackIDEs + = art::ServiceHandle()->HitToTrackIDEs(hit); + float maxe = 0; + float bestfrac = 0; + float bestid = 0; + float ne = 0; + for(size_t i = 0; i < trackIDEs.size(); ++i){ + ne += (float)trackIDEs[i].numElectrons; + if( trackIDEs[i].energy > maxe ) { + maxe = trackIDEs[i].energy; + bestfrac = trackIDEs[i].energyFrac; + bestid = trackIDEs[i].trackID; + } + } + // Save the results + truthid = bestid; + truthidEnergyFrac = bestfrac; + energy = maxe; + numElectrons = ne; + } + + + //================================================================== + // Returns list of all G4 track IDs associated with a hit + std::set HitTruthIds( art::Ptr const& hit){ + std::set ids; + art::ServiceHandle bt_serv; + std::vector trackIDEs = bt_serv->HitToTrackIDEs(hit); + for(size_t i = 0; i < trackIDEs.size(); ++i) ids.insert(trackIDEs[i].trackID); + return ids; + } + */ + + + //===================================================================== + // Get MCTruth associated with TrackID using a try bracket to avoid + // fatal exceptions (return false if no match or exception thrown) + /* + bool G4IdToMCTruth( int const trkID, art::Ptr& mctruth ) + { + art::ServiceHandle pi_serv; + bool matchFound = false; + try { + mctruth = pi_serv->TrackIdToMCTruth_P(trkID); + matchFound = true; + } catch(...) { + std::cout<<"Exception thrown matching TrackID "< geom; + double x = pos.X(); + double y = pos.Y(); + double z = pos.Z(); + auto const& tpcid = geom->FindTPCAtPosition(geo::Point_t{x,y,z}); + if( tpcid.TPC == geo::TPCID::InvalidID ) return -9; + auto const& tpc = geom->TPC(tpcid); + double dx = std::min(x-tpc.MinX(),tpc.MaxX()-x); + double dy = std::min(y-tpc.MinY(),tpc.MaxY()-y); + double dz = std::min(z-tpc.MinZ(),tpc.MaxZ()-z); + return std::min( std::min(dx,dy), dz ); + } + + //=========================================================================== + // Given a line with endpoints L1,L2, return shortest distance betweene the + // line and point 'P' + double DistToLine(TVector3& L1, TVector3& L2, TVector3& p){ + + // general vector formulation: + // a = point on a line + // n = unit vector pointing along line + // --> d = norm[ (p-a) - ((p-a) dot n) * n ] + // In our case, 'a' = L1 + //TVector3 a = L1; + TVector3 n = (L2-L1).Unit(); + TVector3 b = (p-L1); + double projLen = b.Dot(n); + double d = -1; + /* + if ( projLen < 0 ) d = (p-L1).Mag(); + else if ( projLen > (L2-L1).Mag() ) d = (p-L2).Mag(); + else d = (b-projLen*n).Mag(); + */ + if( projLen > 0 && projLen < (L2-L1).Mag() ) { + d = (b-projLen*n).Mag(); + } else { + d = std::min( (p-L1).Mag(), (p-L2).Mag() ); + } + + return d; + } + + double DistToLine2D(TVector2& L1, TVector2& L2, TVector2& p){ + TVector3 newL1(L1.X(), L1.Y(), 0); + TVector3 newL2(L2.X(), L2.Y(), 0); + TVector3 newp(p.X(), p.Y(), 0); + return DistToLine(newL1,newL2,newp); + } + + + //=========================================================================== + bool IsInActiveVolume(geo::Point_t const& p){ + geo::TPCGeo const* TPC = art::ServiceHandle()->PositionToTPCptr(p); + return TPC? TPC->ActiveBoundingBox().ContainsPosition(p): false; + } + + + //========================================================================== + void NormalizeHist(TH1D* h){ + if( h->GetEntries() > 0 ) { + h->Scale(1./h->Integral()); + h->SetBit(TH1::kIsAverage); + h->SetOption("HIST"); + } + } + + + float FindMedian(std::vector& vec){ + if( !vec.size() ) return -9; + size_t n = vec.size() / 2; + std::nth_element(vec.begin(),vec.begin()+n,vec.end()); + if( n % 2 != 0 ) { // odd number of elements + return vec[n]; + }else{ + float a = vec[n]; + std::nth_element(vec.begin(),vec.begin()+n-1,vec.end()); + return (a + vec[n-1]) / 2.0; + } + } + + float FindMean(std::vector& vec){ + float sum = 0; + for(auto& v : vec ) sum += v; + return (vec.size()>0) ? sum/vec.size() : 0; + } + + + /* + //=========================================================================== + float ConvertTicksToX(float peakTime, int plane, int tpc, int cryo) { + auto const* detProp = lar::providerFrom(); + auto const* detClock = lar::providerFrom(); + double Efield = detProp->Efield(0); + double Temp = detProp->Temperature(); + // The drift velocity "Fudge factor"... need to look into this more! + //double fudgeFact = 9.832658e-1; + double driftVel = detProp->DriftVelocity(Efield,Temp)*fudgeFact; + double drift = (peakTime - detProp->GetXTicksOffset(plane,tpc,cryo))*detClock->TPCClock().TickPeriod(); + double X = drift * driftVel; + return X; + } + */ + +} diff --git a/sbnobj/SBND/Blip/BlipUtils.h b/sbnobj/SBND/Blip/BlipUtils.h new file mode 100644 index 00000000..bcd216b3 --- /dev/null +++ b/sbnobj/SBND/Blip/BlipUtils.h @@ -0,0 +1,92 @@ +/////////////////////////////////////////////////////// +// BlipUtils.h +// +// Helper functions and algs. Based on RecoUtils and +// functions in AnalysisTree. +// +////////////////////////////////////////////////////// +#ifndef BLIPUTIL_H_SEEN +#define BLIPUTIL_H_SEEN + +// framework +#include "canvas/Persistency/Common/Ptr.h" +#include "canvas/Persistency/Common/PtrVector.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" + +// LArSoft +#include "lardata/DetectorInfoServices/DetectorClocksService.h" +#include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "nusimdata/SimulationBase/MCParticle.h" +#include "nusimdata/SimulationBase/MCTruth.h" +#include "lardataobj/Simulation/SimEnergyDeposit.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/Track.h" +#include "larsim/MCCheater/BackTrackerService.h" +#include "lardataobj/AnalysisBase/BackTrackerMatchingData.h" +#include "larsim/MCCheater/ParticleInventoryService.h" +#include "larcore/CoreUtils/ServiceUtil.h" +#include "larcore/Geometry/Geometry.h" +#include "larcore/Geometry/WireReadout.h" +//#include "larcorealg/Geometry/GeometryCore.h" +//#include "larcorealg/Geometry/WireReadoutGeom.h" + +// c++ +#include +#include + +#include "sbndcode/BlipRecoSBND/Utils/DataTypes.h" +#include "TH1D.h" + + +typedef std::vector> SEDVec_t; + +geo::View_t kViews[3]={geo::kU, geo::kV, geo::kW}; + +namespace BlipUtils { + + //################################################### + // Functions related to blip reconstruction + //################################################### + //void InitializeDetProps(); + void FillParticleInfo(simb::MCParticle const&, blip::ParticleInfo&, SEDVec_t&, int plane=2); + //void CalcPartDrift(blip::ParticleInfo&, int); + //void CalcTotalDep(float&,int&,float&, SEDVec_t&); + void MakeTrueBlips(std::vector&, std::vector&); + void GrowTrueBlip(blip::ParticleInfo&, blip::TrueBlip&); + void MergeTrueBlips(std::vector&, float); + void GrowHitClust(blip::HitInfo const&, blip::HitClust&); + bool DoHitsOverlap(art::Ptr const&, art::Ptr const&); + bool DoHitClustsOverlap(blip::HitClust const&, blip::HitClust const&); + bool DoHitClustsOverlap(blip::HitClust const&,float,float); + float CalcHitClustsOverlap(blip::HitClust const&, blip::HitClust const&); + float CalcOverlap(float const&, float const&, float const&, float const&); + bool DoHitClustsMatch(blip::HitClust const&, blip::HitClust const&,float); + blip::HitClust MakeHitClust(std::vector const&); + //blip::Blip MakeBlip(std::vector const&); + blip::Blip MakeBlip(std::vector const&, + detinfo::DetectorPropertiesData const&, + detinfo::DetectorClocksData const&); + + int DriftDirX(geo::TPCGeo const&); + + + //################################################### + // General functions + //################################################### + //void HitTruth(art::Ptr const&, int&, float&, float&, float&); + //si_t HitTruthIds( art::Ptr const&); + //bool G4IdToMCTruth( int const, art::Ptr&); + double PathLength(const simb::MCParticle&, TVector3&, TVector3&); + double PathLength(const simb::MCParticle&); + bool IsAncestorOf(int, int, bool); + double DistToBoundary(const recob::Track::Point_t&); + double DistToLine(TVector3&, TVector3&, TVector3&); + double DistToLine2D(TVector2&, TVector2&, TVector2&); + bool IsInActiveVolume(geo::Point_t const&); + void NormalizeHist(TH1D*); + float FindMedian(std::vector&); + float FindMean(std::vector&); + +} + +#endif diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt new file mode 100644 index 00000000..fc847914 --- /dev/null +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -0,0 +1,45 @@ +cet_make_library( + SOURCE + BlipUtils.cc + LIBRARY_NAME + BlipUtils + LIBRARIES + PUBLIC + larcorealg::Geometry + larcore::Geometry_Geometry_service + larsim::Simulation + lardataobj::Simulation + larsim::MCCheater_BackTrackerService_service + larsim::MCCheater_ParticleInventoryService_service + lardata::Utilities + larevt::Filters + lardataobj::RawData + larreco::Calorimetry + lardataobj::RecoBase + lardata::RecoObjects + larpandora::LArPandoraInterface + nusimdata::SimulationBase + cetlib::cetlib + cetlib_except::cetlib_except + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::tfile_support + art_root_io::TFileService_service + art::Persistency_Common + canvas::canvas + art::Persistency_Provenance + art::Utilities + messagefacility::MF_MessageLogger + ROOT::Core + fhiclcpp::fhiclcpp + ROOT::Geom + ROOT::XMLIO + ROOT::Gdml + sbnobj::Common_CRT +) + +art_dictionary(DICTIONARY_LIBRARIES BlipUtils) +install_headers() +install_source() +install_fhicl() diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/DataTypes.h similarity index 99% rename from sbnobj/SBND/Blip/BlipDataTypes.h rename to sbnobj/SBND/Blip/DataTypes.h index 79d7f8dc..023eb99c 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/DataTypes.h @@ -1,3 +1,4 @@ + #include "lardataobj/RecoBase/Hit.h" #include "nusimdata/SimulationBase/MCParticle.h" @@ -169,3 +170,5 @@ namespace blip { }; } + + diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h new file mode 100644 index 00000000..403af300 --- /dev/null +++ b/sbnobj/SBND/Blip/classes.h @@ -0,0 +1,36 @@ +// +// Build a dictionary. +// +// $Id: classes.h,v 1.8 2010/04/12 18:12:28 Exp $ +// $Author: $ +// $Date: 2010/04/12 18:12:28 $ +// +// Original author Rob Kutschke, modified by wes +// + +#include "canvas/Persistency/Common/Wrapper.h" + +// data-products +// lardataobj +//#include "lardata/Utilities/AssociationUtil.h" +#include "canvas/Persistency/Common/Assns.h" +#include "lardataobj/RecoBase/PFParticle.h" +#include "lardataobj/RecoBase/Hit.h" +#include "nusimdata/SimulationBase/MCTruth.h" +#include "sbndcode/BlipRecoSBND/Utils/DataTypes.h" +#include "lardataobj/RecoBase/SpacePoint.h" + +// +// Only include objects that we would like to be able to put into the event. +// Do not include the objects they contain internally. +// + +template class art::Assns; +template class art::Wrapper >; +template class std::vector; +template class art::Wrapper >; +template class std::map; +template class art::Assns; +template class art::Wrapper >; +template class art::Assns; +template class art::Wrapper >; diff --git a/sbnobj/SBND/Blip/classes_def.xml b/sbnobj/SBND/Blip/classes_def.xml new file mode 100644 index 00000000..c5492251 --- /dev/null +++ b/sbnobj/SBND/Blip/classes_def.xml @@ -0,0 +1,14 @@ + + + + + + + + + + + + + + diff --git a/sbnobj/SBND/CMakeLists.txt b/sbnobj/SBND/CMakeLists.txt index a7c6baaf..7eeaf198 100644 --- a/sbnobj/SBND/CMakeLists.txt +++ b/sbnobj/SBND/CMakeLists.txt @@ -3,3 +3,4 @@ add_subdirectory(Commissioning) add_subdirectory(Trigger) add_subdirectory(ToF) add_subdirectory(Timing) +add_subdirectory(Blip) From 9751aa3ef60a30161d85577f6432cb893f4bc8cd Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 10:05:55 -0600 Subject: [PATCH 03/41] Moving blip utils --- sbnobj/SBND/Blip/BlipUtils.h | 2 +- sbnobj/SBND/Blip/classes.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipUtils.h b/sbnobj/SBND/Blip/BlipUtils.h index bcd216b3..a2c77d58 100644 --- a/sbnobj/SBND/Blip/BlipUtils.h +++ b/sbnobj/SBND/Blip/BlipUtils.h @@ -34,7 +34,7 @@ #include #include -#include "sbndcode/BlipRecoSBND/Utils/DataTypes.h" +#include "sbnobj/SBND/Blips//DataTypes.h" #include "TH1D.h" diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index 403af300..2f5ed8c2 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -17,7 +17,7 @@ #include "lardataobj/RecoBase/PFParticle.h" #include "lardataobj/RecoBase/Hit.h" #include "nusimdata/SimulationBase/MCTruth.h" -#include "sbndcode/BlipRecoSBND/Utils/DataTypes.h" +#include "sbnobj/SBND/Blips/DataTypes.h" #include "lardataobj/RecoBase/SpacePoint.h" // From 1387b2a030bb849e612b301116319c1ca3f90921 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 10:10:21 -0600 Subject: [PATCH 04/41] trying to remove depends --- sbnobj/SBND/Blip/CMakeLists.txt | 2 -- 1 file changed, 2 deletions(-) diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index fc847914..fc05f0bd 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -7,8 +7,6 @@ cet_make_library( PUBLIC larcorealg::Geometry larcore::Geometry_Geometry_service - larsim::Simulation - lardataobj::Simulation larsim::MCCheater_BackTrackerService_service larsim::MCCheater_ParticleInventoryService_service lardata::Utilities From f5c315a4e5952ebea0d313be16e1256a8ea433e2 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 10:12:53 -0600 Subject: [PATCH 05/41] trying to remove depends --- sbnobj/SBND/Blip/CMakeLists.txt | 2 -- 1 file changed, 2 deletions(-) diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index fc05f0bd..ff8f2e55 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -7,8 +7,6 @@ cet_make_library( PUBLIC larcorealg::Geometry larcore::Geometry_Geometry_service - larsim::MCCheater_BackTrackerService_service - larsim::MCCheater_ParticleInventoryService_service lardata::Utilities larevt::Filters lardataobj::RawData From 2ce06f9021b145e5a81f6aa6070828d415627b04 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 10:16:06 -0600 Subject: [PATCH 06/41] Im not sure how to link the dictionaries here anymore --- sbnobj/SBND/Blip/BlipUtils.cc | 749 -------------------------------- sbnobj/SBND/Blip/BlipUtils.h | 92 ---- sbnobj/SBND/Blip/CMakeLists.txt | 41 -- 3 files changed, 882 deletions(-) delete mode 100644 sbnobj/SBND/Blip/BlipUtils.cc delete mode 100644 sbnobj/SBND/Blip/BlipUtils.h delete mode 100644 sbnobj/SBND/Blip/CMakeLists.txt diff --git a/sbnobj/SBND/Blip/BlipUtils.cc b/sbnobj/SBND/Blip/BlipUtils.cc deleted file mode 100644 index 1b1e674f..00000000 --- a/sbnobj/SBND/Blip/BlipUtils.cc +++ /dev/null @@ -1,749 +0,0 @@ -#include "BlipUtils.h" - -namespace BlipUtils { - - //============================================================================ - // Find total visible energy deposited in the LAr, and number of electrons deposited - // and drifted to the anode. - /* - void CalcTotalDep(float& energy, int& ne_dep, float& ne_anode, SEDVec_t& sedvec){ - - // energy and electrons deposited - energy = 0; - ne_dep = 0; - for(auto& sed : sedvec ) { - energy += sed->Energy(); - ne_dep += sed->NumElectrons(); - } - - // electrons drifted to collection plane wires - art::ServiceHandle geom; - ne_anode = 0; - for(auto const &chan : art::ServiceHandle()->SimChannels()) { - if( geom->View(chan->Channel()) != geo::kW ) continue; - for(auto const &tdcide : chan->TDCIDEMap() ) { - for(const auto& ide : tdcide.second) ne_anode += ide.numElectrons; - } - } - } - */ - - - //=========================================================================== - // Provided a MCParticle, calculate everything we'll need for later calculations - // and save into ParticleInfo object - void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo, SEDVec_t& sedvec, int caloPlane){ - - // Get important info and do conversions - pinfo.particle = part; - pinfo.trackId = part.TrackId(); - pinfo.isPrimary = (int)(part.Process() == "primary"); - pinfo.mass = /*GeV->MeV*/1e3 * part.Mass(); - pinfo.E = /*GeV->MeV*/1e3 * part.E(); - pinfo.endE = /*GeV->MeV*/1e3 * part.EndE(); - pinfo.KE = /*GeV->MeV*/1e3 * (part.E()-part.Mass()); - pinfo.endKE = /*GeV->MeV*/1e3 * (part.EndE()-part.Mass()); - pinfo.P = /*GeV->MeV*/1e3 * part.Momentum().Vect().Mag(); - pinfo.Px = /*GeV->MeV*/1e3 * part.Px(); - pinfo.Py = /*GeV->MeV*/1e3 * part.Py(); - pinfo.Pz = /*GeV->MeV*/1e3 * part.Pz(); - pinfo.time = /*ns ->mus*/1e-3 * part.T(); - pinfo.endtime = /*ns ->mus*/1e-3 * part.EndT(); - pinfo.numTrajPts = part.NumberTrajectoryPoints(); - - // Pathlength (in AV) and start/end point - pinfo.pathLength = PathLength( part, pinfo.startPoint, pinfo.endPoint); - - // Central position of trajectory - pinfo.position = 0.5*(pinfo.startPoint+pinfo.endPoint); - - // Energy/charge deposited by this particle, found using SimEnergyDeposits - pinfo.depEnergy = 0; - pinfo.depElectrons = 0; - for(auto& sed : sedvec ) { - if( sed->TrackID() == part.TrackId() ) { - pinfo.depEnergy += sed->Energy(); - pinfo.depElectrons += sed->NumElectrons(); - } - } - - return; - - } - - //=================================================================== - // Provided a vector of all particle information for event, fill a - // vector of true blips - void MakeTrueBlips( std::vector& pinfo, std::vector& trueblips ) { - - art::ServiceHandle geom; - auto const& detProp = art::ServiceHandle()->DataForJob(); - auto const& clockData = art::ServiceHandle()->DataForJob(); - - for(size_t i=0; i0) ? tb.Time/clockData.TPCClock().TickPeriod() : 0; - auto point = geo::Point_t{tb.Position.X(),tb.Position.Y(),tb.Position.Z()}; - auto const& tpcID = geom->FindTPCAtPosition(point); - auto const& planeID = geo::PlaneID{tpcID, 0}; - float tick_calc = (float)detProp.ConvertXToTicks(tb.Position.X(),planeID); - tb.DriftTime = tick_calc*clockData.TPCClock().TickPeriod() + clockData.TriggerOffsetTPC(); - - tb.ID = trueblips.size(); - trueblips.push_back(tb); - - } - - } - - - //==================================================================== - void GrowTrueBlip( blip::ParticleInfo& pinfo, blip::TrueBlip& tblip ) { - - simb::MCParticle& part = pinfo.particle; - - // Skip neutrons, photons - if( part.PdgCode() == 2112 || part.PdgCode() == 22 ) return; - - // Check that path length isn't zero - if( !pinfo.pathLength ) return; - - // If this is a new blip, initialize - if( !tblip.G4ChargeMap.size() ) { - tblip.Position = pinfo.position; - tblip.Time = pinfo.time; - - // .. otherwise, check that the new particle - // creation time is comparable to existing blip. - // then calculate new energy-weighted position. - } else if ( fabs(tblip.Time-pinfo.time) < 3 ) { - float totE = tblip.Energy + pinfo.depEnergy; - float w1 = tblip.Energy/totE; - float w2 = pinfo.depEnergy/totE; - tblip.Position = w1*tblip.Position + w2*pinfo.position; - tblip.Time = w1*tblip.Time + w2*pinfo.time; - tblip.LeadCharge = pinfo.depElectrons; - // ... if the particle isn't a match, show's over - } else { - return; - } - - tblip.Energy += pinfo.depEnergy; - tblip.DepElectrons+= pinfo.depElectrons; - tblip.NumElectrons+= std::max(0.,pinfo.numElectrons); - - tblip.G4ChargeMap[part.TrackId()] += pinfo.depElectrons; - tblip.G4PDGMap[part.TrackId()] = part.PdgCode(); - if(pinfo.depElectrons > tblip.LeadCharge ) { - tblip.LeadCharge = pinfo.depElectrons; - tblip.LeadG4Index = pinfo.index; - tblip.LeadG4ID = part.TrackId(); - tblip.LeadG4PDG = part.PdgCode(); - } - } - - - //==================================================================== - // Merge blips that are close - void MergeTrueBlips(std::vector& vtb, float dmin){ - if( dmin <= 0 ) return; - std::vector vtb_merged; - std::vector isGrouped(vtb.size(),false); - - for(size_t i=0; i 5 ) continue; - float d = (blip_i.Position-blip_j.Position).Mag(); - if( d < dmin ) { - isGrouped.at(j) = true; - //float totE = blip_i.Energy + blip_j.Energy; - float totQ = blip_i.DepElectrons + blip_j.DepElectrons; - float w1 = blip_i.DepElectrons/totQ; - float w2 = blip_j.DepElectrons/totQ; - blip_i.Energy += blip_j.Energy; - blip_i.Position = w1*blip_i.Position + w2*blip_j.Position; - blip_i.DriftTime = w1*blip_i.DriftTime+ w2*blip_j.DriftTime; - blip_i.Time = w1*blip_i.Time + w2*blip_j.Time; - blip_i.DepElectrons += blip_j.DepElectrons; - if( blip_j.NumElectrons ) blip_i.NumElectrons += blip_j.NumElectrons; - - blip_i.G4ChargeMap.insert(blip_j.G4ChargeMap.begin(), blip_j.G4ChargeMap.end()); - blip_i.G4PDGMap.insert(blip_j.G4PDGMap.begin(), blip_j.G4PDGMap.end()); - - if( blip_j.LeadCharge > blip_i.LeadCharge ) { - blip_i.LeadCharge = blip_j.LeadCharge; - blip_i.LeadG4ID = blip_j.LeadG4ID; - blip_i.LeadG4Index = blip_j.LeadG4Index; - blip_i.LeadG4PDG = blip_j.LeadG4PDG; - } - }//d < dmin - }//loop over blip_j - blip_i.ID = vtb_merged.size(); - vtb_merged.push_back(blip_i); - } - vtb.clear(); - vtb = vtb_merged; - } - - - //================================================================= - blip::HitClust MakeHitClust(std::vector const& hitinfoVec){ - - blip::HitClust hc; - if( hitinfoVec.size() ) { - int tpc = hitinfoVec[0].tpc; - int plane = hitinfoVec[0].plane; - int cryo = hitinfoVec[0].cryo; - - // check that all hits are on same plane; - for(auto& h : hitinfoVec ) { - if( h.cryo != cryo ) return hc; - if( h.tpc != tpc ) return hc; - if( h.plane != plane ) return hc; - } - - // initialize values - hc.Cryostat = cryo; - hc.TPC = tpc; - hc.Plane = plane; - hc.ADCs = 0; - hc.Charge = 0; - hc.SigmaCharge = 0; - hc.Amplitude = 0; - hc.NPulseTrainHits = 0; - float startTime = 9e9; - float endTime = -9e9; - float weightedTick = 0; - float weightedTime = 0; - float weightedGOF = 0; - //float weightedRatio = 0; - float qGOF = 0; - - // store hit times, charges, and RMS - std::vector tvec; - std::vector qvec; - std::vector dqvec; - std::vector rmsvec; - - // grow our hit cluster! - for(auto& hitinfo : hitinfoVec ) { - if( hc.HitIDs.find(hitinfo.hitid) != hc.HitIDs.end() ) continue; - hc.HitIDs .insert(hitinfo.hitid); - hc.Wires .insert(hitinfo.wire); - hc.Chans .insert(hitinfo.chan); - float q = (hitinfo.charge > 0)? hitinfo.charge : 0; - float integral = hitinfo.integralADC; - float sigma = hitinfo.sigmaintegral; - float dq = (integral != 0 && sigma>0)? (sigma/integral)*q : 0; - hc.Charge += q; - hc.ADCs += hitinfo.integralADC; - hc.Amplitude = std::max(hc.Amplitude, hitinfo.amp ); - weightedTick += q*hitinfo.peakTime; - weightedTime += q*hitinfo.driftTime; - startTime = std::min(startTime, hitinfo.driftTime-hitinfo.rms); - endTime = std::max(endTime, hitinfo.driftTime+hitinfo.rms); - tvec .push_back(hitinfo.driftTime); - qvec .push_back(q); - dqvec .push_back(dq); - rmsvec .push_back(hitinfo.rms); - if( hitinfo.g4trkid >= 0 ) hc.G4IDs.insert(hitinfo.g4trkid); - if( hitinfo.gof < 0 ) { - hc.NPulseTrainHits++; - } else { - weightedGOF += q*hitinfo.gof; - qGOF += q; - } - }//endloop over hits - - // mean goodness of fit - if( qGOF ) hc.GoodnessOfFit = weightedGOF/qGOF; - - // calculate other quantities - hc.NHits = hc.HitIDs.size(); - hc.NWires = hc.Wires.size(); - hc.CenterWire =(*hc.Wires.begin()+*hc.Wires.rbegin())/2.; - hc.CenterChan =(*hc.Chans.begin()+*hc.Chans.rbegin())/2.; - hc.StartWire = *hc.Wires.begin(); - hc.EndWire = *hc.Wires.rbegin(); - hc.StartTime = startTime; - hc.EndTime = endTime; - hc.Timespan = hc.EndTime - hc.StartTime; - hc.Time = weightedTime / hc.Charge; - hc.TimeTick = weightedTick / hc.Charge; - - // overall cluster RMS and uncertainty in charge - float sig_sumSq = 0; - float dt_sumSq = 0; - float dq = 0; - for(size_t i=0; i 0 hits - - // mark the cluster as valid and ship it out - hc.isValid = true; - return hc; - } - - - //================================================================= - blip::Blip MakeBlip( std::vector const& hcs, - detinfo::DetectorPropertiesData const& detProp, - detinfo::DetectorClocksData const& clockData ){ - - art::ServiceHandle wireReadoutGeom; - blip::Blip newblip; - - // ------------------------------------------------ - // Must be 1-3 clusts (no more, no less!) - if( hcs.size() > 3 || hcs.size() < 1 ) return newblip; - - // ------------------------------------------------ - // All hits must be in same TPC, and no 2 planes can be repeated - std::set planeIDs; - for(size_t i=0; i()->DataForJob(); - //auto const& clockData = art::ServiceHandle()->DataForJob(); - float driftVelocity = detProp.DriftVelocity(detProp.Efield(0),detProp.Temperature()); - float tick_to_cm = clockData.TPCClock().TickPeriod() * driftVelocity; - - - newblip.Cryostat = hcs[0].Cryostat; - newblip.TPC = hcs[0].TPC; - newblip.NPlanes = planeIDs.size(); - - // ------------------------------------------------ - /// Look for valid wire intersections between - // central-most hits in each cluster - std::vector wirex; - for(size_t i=0; iGet().Plane(geo::PlaneID{(unsigned int)hcs[i].Cryostat, (unsigned int)hcs[i].TPC, (unsigned int)hcs[i].Plane}); - double wirepitch = planegeo.WirePitch(); - // use view with the maximal wire extent to calculate transverse (YZ) length - if( hcs[i].NWires > newblip.MaxWireSpan ) { - newblip.MaxWireSpan = hcs[i].NWires; - newblip.dYZ = hcs[i].NWires * wirepitch; - } - - for(size_t j=i+1; jsecond.Y()); - intsec_p.SetZ(hcs[i].IntersectLocations.find(hcs[j].ID)->second.Z()); - } else { - std::vector i_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcs[i].CenterChan); - std::vector j_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcs[j].CenterChan); - - match3d = wireReadoutGeom->Get().WireIDsIntersect(i_wireids.at(0), j_wireids.at(0), intsec_p); - } - - if( match3d ) { - TVector3 a(0., intsec_p.Y(), intsec_p.Z()); - wirex.push_back(a); - newblip.clusters[pli] = hcs[i]; - newblip.clusters[plj] = hcs[j]; - } - } - } - - // Require some number of intersection points. - if( !wirex.size() ) return newblip; - - // Loop over the intersection points and calculate average position in - // YZ-plane, as well as the mean difference between intersection points. - newblip.Position.SetXYZ(0,0,0); - if( wirex.size() == 1 ) { - newblip.Position= wirex[0]; - } else { - newblip.SigmaYZ = 0; - double fact = 1./wirex.size(); - for(auto& v : wirex ) newblip.Position += v * fact; - for(auto& v : wirex ) newblip.SigmaYZ += (v-newblip.Position).Mag() * fact; - // Ensure that difference between intersection points is - // consistent with the maximal wire extent - if( newblip.SigmaYZ > std::max(1.,0.5*newblip.dYZ) ) return newblip; - } - - // Calculate mean drift time and X-position - // (note that the 'time' of each of the hit clusters - // have already been corrected for plane-to-plane offsets) - newblip.TimeTick= 0; - newblip.Time = 0; - newblip.dX = 0; - float vsize = (float)hcs.size(); - for(auto hc : hcs ) { - newblip.TimeTick += hc.TimeTick / vsize; - newblip.Time += hc.Time / vsize; - newblip.dX = std::max((float)(hc.EndTime-hc.StartTime)*tick_to_cm, newblip.dX); - } - - //auto const& tpcID = geo::TPCID(geo::CryostatID(newblip.Cryostat),newblip.TPC); - //auto const& planeID = art::ServiceHandle()->GetBeginPlaneID(tpcID); - //newblip.Position .SetX(detProp.ConvertTicksToX(newblip.TimeTick, planeID)); - - // convert ticks to X - auto const& cryostat= art::ServiceHandle()->Cryostat(geo::CryostatID(newblip.Cryostat)); - auto const& tpcgeom = cryostat.TPC(newblip.TPC); - auto tpcID = tpcgeom.ID(); - auto const& planegeo = wireReadoutGeom->Get().Plane(geo::PlaneID{tpcID, 0}); - auto const xyz = planegeo.GetCenter(); - int dirx = DriftDirX(tpcgeom); - - newblip.Position.SetX( xyz.X() + dirx * tick_to_cm * newblip.Time ); - - // this should ALREADY be accounted for at the hit-processing level in BlipRecoAlg, - // through the use of GetXTicksOffset... - //float offset_ticks = clockData.TriggerOffsetTPC() / clockData.TPCClock().TickPeriod(); - //float driftTicks = newblip.Time + clockData.TriggerOffsetTPC(); - //std::cout<<"blip time "< pi_serv; - const sim::ParticleList& plist = pi_serv->ParticleList(); - if( particleID == ancestorID ) return true; - if( particleID < ancestorID ) return false; - if( !plist.HasParticle(ancestorID) ) return false; - while( particleID > ancestorID ) { - - simb::MCParticle p = pi_serv->TrackIdToParticle(particleID); - if( !plist.HasParticle(p.Mother() ) ) { return false; } - - simb::MCParticle pM = pi_serv->TrackIdToParticle(p.Mother()); - if ( pM.TrackId() == ancestorID ) { return true; } - else if ( breakAtPhots == true && pM.PdgCode() == 22 ) { return false; } - else if ( pM.Process() == "primary" || pM.TrackId() == 1 ) { return false; } - else if ( pM.Mother() == 0 ) { return false; } - else { particleID = pM.TrackId(); } - } - - return false; - } - - //=================================================================== - int DriftDirX(geo::TPCGeo const& tpcgeom) { - return ((tpcgeom.DriftSign() == geo::DriftSign::Negative) ? +1.0 : -1.0); - } - - //==================================================================== - bool DoHitsOverlap(art::Ptr const& hit1, art::Ptr const& hit2){ - if( hit1->WireID() != hit2->WireID() ) return false; - float t1 = hit1->PeakTime(); - float t2 = hit2->PeakTime(); - float sig = std::max(hit1->RMS(),hit2->RMS()); - if( fabs(t1-t2) < 1.0*sig ) return true; - else return false; - } - - - //==================================================================== - bool DoHitClustsOverlap(blip::HitClust const& hc1, blip::HitClust const& hc2){ - - // only match across different wires in same TPC - if( hc1.TPC != hc2.TPC ) return false; - - if( hc1.StartTime <= hc2.EndTime - && hc2.StartTime <= hc1.EndTime ) return true; - else return false; - } - bool DoHitClustsOverlap(blip::HitClust const& hc1, float t1, float t2 ){ - blip::HitClust hc2; - hc2.TPC = hc1.TPC; - hc2.StartTime = t1; - hc2.EndTime = t2; - return DoHitClustsOverlap(hc1,hc2); - } - - //==================================================================== - // Calculates the level of time overlap between two clusters - float CalcHitClustsOverlap(blip::HitClust const& hc1, blip::HitClust const& hc2){ - return CalcOverlap(hc1.StartTime,hc1.EndTime,hc2.StartTime,hc2.EndTime); - } - - float CalcOverlap(const float& x1, const float& x2, const float& y1, const float& y2){ - float full_range = std::max(x2,y2) - std::min(x1,y1); - float sum = (x2-x1) + (y2-y1); - float overlap = std::max(float(0), sum-full_range); - if( overlap > 0 ) return 2. * overlap / sum; - else return -1; - } - - //==================================================================== - bool DoHitClustsMatch(blip::HitClust const& hc1, blip::HitClust const& hc2, float minDiffTicks = 2){ - if( fabs(hc1.Time-hc2.Time) < minDiffTicks ) return true; - else return false; - } - - //==================================================================== - // This function calculates the leading MCParticle ID contributing to a hit and the - // fraction of that hit's energy coming from that particle. - /* - void HitTruth(art::Ptr const& hit, int& truthid, float& truthidEnergyFrac, float& energy,float& numElectrons){ - // Get associated sim::TrackIDEs for this hit - std::vector trackIDEs - = art::ServiceHandle()->HitToTrackIDEs(hit); - float maxe = 0; - float bestfrac = 0; - float bestid = 0; - float ne = 0; - for(size_t i = 0; i < trackIDEs.size(); ++i){ - ne += (float)trackIDEs[i].numElectrons; - if( trackIDEs[i].energy > maxe ) { - maxe = trackIDEs[i].energy; - bestfrac = trackIDEs[i].energyFrac; - bestid = trackIDEs[i].trackID; - } - } - // Save the results - truthid = bestid; - truthidEnergyFrac = bestfrac; - energy = maxe; - numElectrons = ne; - } - - - //================================================================== - // Returns list of all G4 track IDs associated with a hit - std::set HitTruthIds( art::Ptr const& hit){ - std::set ids; - art::ServiceHandle bt_serv; - std::vector trackIDEs = bt_serv->HitToTrackIDEs(hit); - for(size_t i = 0; i < trackIDEs.size(); ++i) ids.insert(trackIDEs[i].trackID); - return ids; - } - */ - - - //===================================================================== - // Get MCTruth associated with TrackID using a try bracket to avoid - // fatal exceptions (return false if no match or exception thrown) - /* - bool G4IdToMCTruth( int const trkID, art::Ptr& mctruth ) - { - art::ServiceHandle pi_serv; - bool matchFound = false; - try { - mctruth = pi_serv->TrackIdToMCTruth_P(trkID); - matchFound = true; - } catch(...) { - std::cout<<"Exception thrown matching TrackID "< geom; - double x = pos.X(); - double y = pos.Y(); - double z = pos.Z(); - auto const& tpcid = geom->FindTPCAtPosition(geo::Point_t{x,y,z}); - if( tpcid.TPC == geo::TPCID::InvalidID ) return -9; - auto const& tpc = geom->TPC(tpcid); - double dx = std::min(x-tpc.MinX(),tpc.MaxX()-x); - double dy = std::min(y-tpc.MinY(),tpc.MaxY()-y); - double dz = std::min(z-tpc.MinZ(),tpc.MaxZ()-z); - return std::min( std::min(dx,dy), dz ); - } - - //=========================================================================== - // Given a line with endpoints L1,L2, return shortest distance betweene the - // line and point 'P' - double DistToLine(TVector3& L1, TVector3& L2, TVector3& p){ - - // general vector formulation: - // a = point on a line - // n = unit vector pointing along line - // --> d = norm[ (p-a) - ((p-a) dot n) * n ] - // In our case, 'a' = L1 - //TVector3 a = L1; - TVector3 n = (L2-L1).Unit(); - TVector3 b = (p-L1); - double projLen = b.Dot(n); - double d = -1; - /* - if ( projLen < 0 ) d = (p-L1).Mag(); - else if ( projLen > (L2-L1).Mag() ) d = (p-L2).Mag(); - else d = (b-projLen*n).Mag(); - */ - if( projLen > 0 && projLen < (L2-L1).Mag() ) { - d = (b-projLen*n).Mag(); - } else { - d = std::min( (p-L1).Mag(), (p-L2).Mag() ); - } - - return d; - } - - double DistToLine2D(TVector2& L1, TVector2& L2, TVector2& p){ - TVector3 newL1(L1.X(), L1.Y(), 0); - TVector3 newL2(L2.X(), L2.Y(), 0); - TVector3 newp(p.X(), p.Y(), 0); - return DistToLine(newL1,newL2,newp); - } - - - //=========================================================================== - bool IsInActiveVolume(geo::Point_t const& p){ - geo::TPCGeo const* TPC = art::ServiceHandle()->PositionToTPCptr(p); - return TPC? TPC->ActiveBoundingBox().ContainsPosition(p): false; - } - - - //========================================================================== - void NormalizeHist(TH1D* h){ - if( h->GetEntries() > 0 ) { - h->Scale(1./h->Integral()); - h->SetBit(TH1::kIsAverage); - h->SetOption("HIST"); - } - } - - - float FindMedian(std::vector& vec){ - if( !vec.size() ) return -9; - size_t n = vec.size() / 2; - std::nth_element(vec.begin(),vec.begin()+n,vec.end()); - if( n % 2 != 0 ) { // odd number of elements - return vec[n]; - }else{ - float a = vec[n]; - std::nth_element(vec.begin(),vec.begin()+n-1,vec.end()); - return (a + vec[n-1]) / 2.0; - } - } - - float FindMean(std::vector& vec){ - float sum = 0; - for(auto& v : vec ) sum += v; - return (vec.size()>0) ? sum/vec.size() : 0; - } - - - /* - //=========================================================================== - float ConvertTicksToX(float peakTime, int plane, int tpc, int cryo) { - auto const* detProp = lar::providerFrom(); - auto const* detClock = lar::providerFrom(); - double Efield = detProp->Efield(0); - double Temp = detProp->Temperature(); - // The drift velocity "Fudge factor"... need to look into this more! - //double fudgeFact = 9.832658e-1; - double driftVel = detProp->DriftVelocity(Efield,Temp)*fudgeFact; - double drift = (peakTime - detProp->GetXTicksOffset(plane,tpc,cryo))*detClock->TPCClock().TickPeriod(); - double X = drift * driftVel; - return X; - } - */ - -} diff --git a/sbnobj/SBND/Blip/BlipUtils.h b/sbnobj/SBND/Blip/BlipUtils.h deleted file mode 100644 index a2c77d58..00000000 --- a/sbnobj/SBND/Blip/BlipUtils.h +++ /dev/null @@ -1,92 +0,0 @@ -/////////////////////////////////////////////////////// -// BlipUtils.h -// -// Helper functions and algs. Based on RecoUtils and -// functions in AnalysisTree. -// -////////////////////////////////////////////////////// -#ifndef BLIPUTIL_H_SEEN -#define BLIPUTIL_H_SEEN - -// framework -#include "canvas/Persistency/Common/Ptr.h" -#include "canvas/Persistency/Common/PtrVector.h" -#include "art/Framework/Services/Registry/ServiceHandle.h" - -// LArSoft -#include "lardata/DetectorInfoServices/DetectorClocksService.h" -#include "lardata/DetectorInfoServices/DetectorPropertiesService.h" -#include "nusimdata/SimulationBase/MCParticle.h" -#include "nusimdata/SimulationBase/MCTruth.h" -#include "lardataobj/Simulation/SimEnergyDeposit.h" -#include "lardataobj/RecoBase/Hit.h" -#include "lardataobj/RecoBase/Track.h" -#include "larsim/MCCheater/BackTrackerService.h" -#include "lardataobj/AnalysisBase/BackTrackerMatchingData.h" -#include "larsim/MCCheater/ParticleInventoryService.h" -#include "larcore/CoreUtils/ServiceUtil.h" -#include "larcore/Geometry/Geometry.h" -#include "larcore/Geometry/WireReadout.h" -//#include "larcorealg/Geometry/GeometryCore.h" -//#include "larcorealg/Geometry/WireReadoutGeom.h" - -// c++ -#include -#include - -#include "sbnobj/SBND/Blips//DataTypes.h" -#include "TH1D.h" - - -typedef std::vector> SEDVec_t; - -geo::View_t kViews[3]={geo::kU, geo::kV, geo::kW}; - -namespace BlipUtils { - - //################################################### - // Functions related to blip reconstruction - //################################################### - //void InitializeDetProps(); - void FillParticleInfo(simb::MCParticle const&, blip::ParticleInfo&, SEDVec_t&, int plane=2); - //void CalcPartDrift(blip::ParticleInfo&, int); - //void CalcTotalDep(float&,int&,float&, SEDVec_t&); - void MakeTrueBlips(std::vector&, std::vector&); - void GrowTrueBlip(blip::ParticleInfo&, blip::TrueBlip&); - void MergeTrueBlips(std::vector&, float); - void GrowHitClust(blip::HitInfo const&, blip::HitClust&); - bool DoHitsOverlap(art::Ptr const&, art::Ptr const&); - bool DoHitClustsOverlap(blip::HitClust const&, blip::HitClust const&); - bool DoHitClustsOverlap(blip::HitClust const&,float,float); - float CalcHitClustsOverlap(blip::HitClust const&, blip::HitClust const&); - float CalcOverlap(float const&, float const&, float const&, float const&); - bool DoHitClustsMatch(blip::HitClust const&, blip::HitClust const&,float); - blip::HitClust MakeHitClust(std::vector const&); - //blip::Blip MakeBlip(std::vector const&); - blip::Blip MakeBlip(std::vector const&, - detinfo::DetectorPropertiesData const&, - detinfo::DetectorClocksData const&); - - int DriftDirX(geo::TPCGeo const&); - - - //################################################### - // General functions - //################################################### - //void HitTruth(art::Ptr const&, int&, float&, float&, float&); - //si_t HitTruthIds( art::Ptr const&); - //bool G4IdToMCTruth( int const, art::Ptr&); - double PathLength(const simb::MCParticle&, TVector3&, TVector3&); - double PathLength(const simb::MCParticle&); - bool IsAncestorOf(int, int, bool); - double DistToBoundary(const recob::Track::Point_t&); - double DistToLine(TVector3&, TVector3&, TVector3&); - double DistToLine2D(TVector2&, TVector2&, TVector2&); - bool IsInActiveVolume(geo::Point_t const&); - void NormalizeHist(TH1D*); - float FindMedian(std::vector&); - float FindMean(std::vector&); - -} - -#endif diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt deleted file mode 100644 index ff8f2e55..00000000 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ /dev/null @@ -1,41 +0,0 @@ -cet_make_library( - SOURCE - BlipUtils.cc - LIBRARY_NAME - BlipUtils - LIBRARIES - PUBLIC - larcorealg::Geometry - larcore::Geometry_Geometry_service - lardata::Utilities - larevt::Filters - lardataobj::RawData - larreco::Calorimetry - lardataobj::RecoBase - lardata::RecoObjects - larpandora::LArPandoraInterface - nusimdata::SimulationBase - cetlib::cetlib - cetlib_except::cetlib_except - art::Framework_Core - art::Framework_Principal - art::Framework_Services_Registry - art_root_io::tfile_support - art_root_io::TFileService_service - art::Persistency_Common - canvas::canvas - art::Persistency_Provenance - art::Utilities - messagefacility::MF_MessageLogger - ROOT::Core - fhiclcpp::fhiclcpp - ROOT::Geom - ROOT::XMLIO - ROOT::Gdml - sbnobj::Common_CRT -) - -art_dictionary(DICTIONARY_LIBRARIES BlipUtils) -install_headers() -install_source() -install_fhicl() From 2d37c3fb1dfe7f546378f26bdc60581d4c11bb27 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 10:19:40 -0600 Subject: [PATCH 07/41] Simplifying includes --- sbnobj/SBND/Blip/CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 sbnobj/SBND/Blip/CMakeLists.txt diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt new file mode 100644 index 00000000..88b49a2e --- /dev/null +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -0,0 +1,3 @@ +install_headers() +install_source() +install_fhicl() From 55d9f70ffd9fac8afda1445c91e9eb9407cfcfe0 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 11:28:42 -0600 Subject: [PATCH 08/41] trying to make a library --- sbnobj/SBND/Blip/{DataTypes.h => BlipDataTypes.h} | 0 sbnobj/SBND/Blip/CMakeLists.txt | 1 + 2 files changed, 1 insertion(+) rename sbnobj/SBND/Blip/{DataTypes.h => BlipDataTypes.h} (100%) diff --git a/sbnobj/SBND/Blip/DataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h similarity index 100% rename from sbnobj/SBND/Blip/DataTypes.h rename to sbnobj/SBND/Blip/BlipDataTypes.h diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index 88b49a2e..82dd8b31 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -1,3 +1,4 @@ +art_dictionary(DICTIONARY_LIBRARIES sbnobj::BlipDataTypes) install_headers() install_source() install_fhicl() From 343d755b539a2f909caa671ca7138664d860688a Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 11:34:12 -0600 Subject: [PATCH 09/41] Still trying to make a library --- sbnobj/SBND/Blip/BlipDataTypes.cc | 1 + sbnobj/SBND/Blip/CMakeLists.txt | 6 ++++++ 2 files changed, 7 insertions(+) create mode 100644 sbnobj/SBND/Blip/BlipDataTypes.cc diff --git a/sbnobj/SBND/Blip/BlipDataTypes.cc b/sbnobj/SBND/Blip/BlipDataTypes.cc new file mode 100644 index 00000000..6006e9f4 --- /dev/null +++ b/sbnobj/SBND/Blip/BlipDataTypes.cc @@ -0,0 +1 @@ +#include "sbnobj/SBND/Blip/BlipDataTypes.cc" diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index 82dd8b31..67a5186c 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -1,3 +1,9 @@ +cet_make_library( + SOURCE + BlipDataTypes.cc + LIBRARIES + ROOT::Core +) art_dictionary(DICTIONARY_LIBRARIES sbnobj::BlipDataTypes) install_headers() install_source() From 99a3c011bcd6ed367756457ea014ae46b41fa3cd Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 6 Nov 2025 12:38:35 -0600 Subject: [PATCH 10/41] sorting compiles --- sbnobj/SBND/Blip/BlipDataTypes.cc | 2 +- sbnobj/SBND/Blip/BlipDataTypes.h | 4 ++-- sbnobj/SBND/Blip/CMakeLists.txt | 7 ++++--- sbnobj/SBND/Blip/classes.h | 2 +- 4 files changed, 8 insertions(+), 7 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.cc b/sbnobj/SBND/Blip/BlipDataTypes.cc index 6006e9f4..69eb2f99 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.cc +++ b/sbnobj/SBND/Blip/BlipDataTypes.cc @@ -1 +1 @@ -#include "sbnobj/SBND/Blip/BlipDataTypes.cc" +#include "sbnobj/SBND/Blip/BlipDataTypes.h" diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 023eb99c..450d7a9f 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -1,9 +1,9 @@ #include "lardataobj/RecoBase/Hit.h" #include "nusimdata/SimulationBase/MCParticle.h" - #include - +#include +#include typedef std::vector vint_t; typedef std::vector vbool_t; typedef std::vector vfloat_t; diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index 67a5186c..60536524 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -1,10 +1,11 @@ -cet_make_library( +cet_make_library( LIBRARY_NAME sbndobj_BlipDataTypes SOURCE BlipDataTypes.cc LIBRARIES - ROOT::Core + lardataobj::RecoBase + nusimdata::SimulationBase ) -art_dictionary(DICTIONARY_LIBRARIES sbnobj::BlipDataTypes) +art_dictionary(DICTIONARY_LIBRARIES sbndobj_BlipDataTypes) install_headers() install_source() install_fhicl() diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index 2f5ed8c2..7508db02 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -17,7 +17,7 @@ #include "lardataobj/RecoBase/PFParticle.h" #include "lardataobj/RecoBase/Hit.h" #include "nusimdata/SimulationBase/MCTruth.h" -#include "sbnobj/SBND/Blips/DataTypes.h" +#include "sbnobj/SBND/Blip/BlipDataTypes.h" #include "lardataobj/RecoBase/SpacePoint.h" // From bcbb2aa46b34886f8b104a40bdf7efc439e2c2c2 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Mon, 10 Nov 2025 15:28:59 -0600 Subject: [PATCH 11/41] compiled version --- sbnobj/SBND/Blip/BlipDataTypes.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 450d7a9f..dc5a5336 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -10,10 +10,10 @@ typedef std::vector vfloat_t; typedef std::set si_t; typedef std::map mif_t; -const int kNplanes = 3; -const int kNTPCs = 2; - namespace blip { + const int kNplanes = 3; + const int kNTPCs = 2; + //################################################### // Data structures From 5bfaf2c6baba60b541b4487a0894f1aa6b210cac Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Thu, 13 Nov 2025 18:47:52 -0600 Subject: [PATCH 12/41] small fix for sbndcode --- sbnobj/SBND/Blip/BlipDataTypes.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index dc5a5336..c59004e7 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -1,4 +1,5 @@ - +#ifndef BLIPDATATYPE +#define BLIPDATATYPE #include "lardataobj/RecoBase/Hit.h" #include "nusimdata/SimulationBase/MCParticle.h" #include @@ -170,5 +171,4 @@ namespace blip { }; } - - +#endif From a535dc9cab54757a67a03f4c4e2358e3e1dac8d7 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Fri, 21 Nov 2025 15:38:10 -0600 Subject: [PATCH 13/41] Small fix for sbndcode addition --- sbnobj/SBND/Blip/classes_def.xml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/sbnobj/SBND/Blip/classes_def.xml b/sbnobj/SBND/Blip/classes_def.xml index c5492251..b6cf251f 100644 --- a/sbnobj/SBND/Blip/classes_def.xml +++ b/sbnobj/SBND/Blip/classes_def.xml @@ -3,6 +3,8 @@ + + From bbd96f789157b1b6daeb261736326a6682c72b0a Mon Sep 17 00:00:00 2001 From: root Date: Sat, 29 Nov 2025 11:24:13 -0600 Subject: [PATCH 14/41] Added lots of doxygen comments. Need to test geo::Point_t still --- sbnobj/SBND/Blip/BlipDataTypes.h | 339 +++++++++++++++++-------------- sbnobj/SBND/Blip/CMakeLists.txt | 1 + 2 files changed, 186 insertions(+), 154 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index c59004e7..94d1e546 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -2,173 +2,204 @@ #define BLIPDATATYPE #include "lardataobj/RecoBase/Hit.h" #include "nusimdata/SimulationBase/MCParticle.h" +#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h" #include #include #include -typedef std::vector vint_t; -typedef std::vector vbool_t; -typedef std::vector vfloat_t; -typedef std::set si_t; -typedef std::map mif_t; namespace blip { - const int kNplanes = 3; - const int kNTPCs = 2; + const int kNplanes = 3; ///< Wire planes in a given tpc in this detector + const int kNTPCs = 2; ///< TPC in this detector - - //################################################### - // Data structures - //################################################### - - struct ParticleInfo { - simb::MCParticle particle; - int trackId = -9; - int index = -9; - int isPrimary = -9; - int numTrajPts = -9; - double depEnergy = -9; - int depElectrons = -9; - double numElectrons = -9; - double mass = -9; - double E = -9; - double endE = -9; - double KE = -9; - double endKE = -9; - double P = -9; - double Px = -9; - double Py = -9; - double Pz = -9; - double pathLength = -9; - double time = -9; - double endtime = -9; - TVector3 startPoint; - TVector3 endPoint; - TVector3 position; +/** Blips are small, plane-matched, energy deposition in liquid argon with sizes similar to wire separation. + * Blips are composed of hit-clusters, which are (time and space) adjacent hits on a single wire plane. + * A single plane, usually the collection, is given a privledged position as the calorimetry plane. + * Every hit-cluster on the calorimetry plane is checked for matches on the other two planes, and when multiple possible pairs are found + * the highest scoring one is selected. + * Score components include cluster charge, cluster time duration, cluster peak time, and wire intersections. + */ + struct Blip { + int ID = -9; ///< Internal index for blip labelling + bool isValid = false; ///< Blip passes basic checks in reco. All blips saved to artRoot file should be valid. + int Cryostat = -9; ///< Which cryostat the blip was reconstructed to. For SBND this should always be 0. + int TPC = -9; ///< Which tpc the blip was reconstructed to. For SBND this may be 0 or 1. + int NPlanes = -9; ///< Number of planes matched to build the blip. Blips must be plane matched so this should be 2+ planes. + int MaxWireSpan = -9; ///< Maximum span of wires on any plane hit-cluster + float TimeTick = -999; ///< Average time of hitclusters making up blips. [tpc tick] + /*! + Hit cluster time is the charge-weighted average of the hit-peak times for hits in the cluster + */ + float Time = -999; ///< Average time of hitclusters making up blips. [us] + /*! + Hit cluster time is the charge-weighted average of the hit-peak times for hits in the cluster + */ + float Charge = -9; ///< Charge on calorimetry plane [e-] + float Energy = -999; ///< Reconstructed energy in the calorimetry plane (const dE/dx, fcl-configurable) [MeV] + float EnergyESTAR = -999; ///< Reconstructed energy in the calorimetry plane (ESTAR method from ArgoNeuT) [MeV] + float EnergyPSTAR = -999; ///< Reconstructed energy in the calorimetry plane (PSTAR method similar with ESTAR method from ArgoNeuT) [MeV] + float ProxTrkDist = -9; ///< 3-D distance to closest track, assuming the blip was concident with event trigger [cm] + /*! + for properly flash matched out-of-time tracks this distance will be wrong! The blips have no such flash matching ability as of yet + */ + int ProxTrkID = -9; ///< index of the of closest track, assuming the blip was concident with event trigger [cm] + bool inCylinder = false; ///< Bool for whether this blip falls in a 45 degree cone relative to any track closer than fcl-set CylinderRadius (15 cm) + /*! + please note the blip X position is unreliable, so these distance and 3-d position derived variables may be incorrect + */ + geo::Point_t Position; ///< 3D position TVector3 //Used to be TVector3 + float SigmaYZ = -9.; ///< Uncertainty in YZ intersect [cm] + float dX = -9; ///< Equivalent length along drift direction [cm] + float dYZ = -9; ///< Approximate length scale in YZ space [cm] + blip::HitClust clusters[kNplanes]; ///< Plane/cluster-specific information. Just contains hit clusters making up this blip! + blip::TrueBlip truth; ///< Truth-matched energy deposition. Just contains the relevant MC truth info to this blip! + double X() { return Position.X(); } + double Y() { return Position.Y(); } + double Z() { return Position.Z(); } + }; - - // True energy depositions - struct TrueBlip { - int ID = -9; // unique blip ID - int Cryostat = -9; // Cryostat ID - int TPC = -9; // TPC ID - float Time = -9; // time of particle interaction - int TimeTick = -9; // time tick - float DriftTime = -9; // drift time [us] - float Energy = 0; // energy dep [MeV] - int DepElectrons = 0; // deposited electrons - int NumElectrons = 0; // electrons reaching wires - int LeadG4ID = -9; // lead G4 track ID - int LeadG4Index = -9; // lead G4 track index - int LeadG4PDG = -9; // lead G4 PDG - float LeadCharge = -9; // lead G4 charge dep - TVector3 Position; // XYZ position - mif_t G4ChargeMap; - mif_t G4PDGMap; + + /** Hit clusters are collections of adjacent hits on a single wire plane. + * They have a fcl-set maximum wire-span (fcl set. Default 10), as well as maximum tick-width (fcl set. Default 50 tick). + * Hit clusters have timestamp and associated wire IDs. + * Within a hit cluster certain statistical summaries of the collection are saved including: + * Total charge, total charge uncertianty, peak hit amplitude, and charge weighted RMS hit spread + */ + struct HitClust { + int ID = -9; ///< Per-plane index for the hit clusters. In SBND we save every collection plane hitcluster but not the induction + bool isValid = false; ///< Bool check that every hit is in the same cryostat, tpc, plane. Should always be true for saved items + int CenterChan = -999; ///< Channel ID of the wire in the geometric center of the hit cluster + int CenterWire = -999; ///< Wire ID of the wire in the geometric center of the hit cluster + bool isTruthMatched = false; ///< is there a trueBlip with the same leadG4ID as one of the G4IDs making up this cluster + bool isMerged = false; ///< Depreciated internal flag for tracking blips before/after merge. + bool isMatched = false; ///< Is this hit cluster plane-matched into a full 3d blip + int DeadWireSep = 99; ///< Separation between the extreme ends of the hitcluster and the nearest dead wire. + /*! + DeadWireSep can be between 0 and 5 and valid. Larger separations are filled in as 99. + */ + int Cryostat = -9; ///< cryostat for this hit cluster + int TPC = -9; ///< TPC for this hit cluster + int Plane = -9; ///< Plane index for this hit cluster + int NHits = -9; ///< Number of hits making up this hit cluster + int NWires = -9; ///< Wire span of the hit cluster + float ADCs = -999; ///< ADC integral sum of hits making up the cluster [ADC-tick] + float Amplitude = -999; ///< Max amplitude of hits making up the hit cluster [ADC] + float Charge = -999; ///< Total charge of hits making up the hit cluster [e-] + /*! + Charge is reconstructed from calo::CalorimetryAlg ( "sbnd_calorimetryalgmc" )-> ElectronsFromADCArea function + Configuration is in sbndcode/sbndcode/LArSoftConfigurations/calorimetry_sbnd.fcl + */ + float SigmaCharge = -999; ///< charge-weighted charge uncertainties for this hit-cluster [e-] + float TimeTick = -999; ///< charge-weighted average hit-peak-times for this hit-cluster [tick] + float Time = -999; ///< charge-weighted average hit-peak-times for this hit-cluster [us] + float StartHitTime = -999; ///< Depreciated + float EndHitTime = -999; ///< Depreciated + float StartTime = -999; ///< Minimum -1 sigma time of a hit in this cluster [us] + float EndTime = -999; ///< Max +1 sigma time of a hit in this cluster [us] + float Timespan = -999; ///< Hit cluster EndTime - StartTime [us] + float RMS = -999; ///< Quadrature estimate of charge waveform timespread accounting for varied hit-drift times and internal hit-RMS [us] + int StartWire = -999; ///< Lowest wireID involved with the hitcluster + int EndWire = -999; ///< Highest wireID involved with the hit cluster + int NPulseTrainHits = -9; ///< Number of hits with a GoodnessOfFit<0 involved in this hit cluster + float GoodnessOfFit = -999; ///< Charge weighted hit-GoodnessofFit param + int BlipID = -9; ///< If this hit cluster ended up in a blip, what is its ID + int EdepID = -9; ///< If this hit cluster is MC-matched what is the trueBlip ID + std::set HitIDs; ///< Index of the recob::hit objects making up this cluster + std::set Wires; ///< Set of geo::wireIDs contributing hits to this cluster + std::set Chans; ///< Set of raw::ChannelID_t contributing hits to this cluster + std::set G4IDs; ///< simb::MCParticle track ID contributing hits to this cluster + std::map IntersectLocations; ///< Internal reconstruction variable for recording where hit-clusters on different planes would overlap //used to be TVector3 }; + /** Extra struct used in blip reco for grabbing key recob::hit information + * Most attributes are grabbed directly from the recob::hit, but some, like charge, involve extra processing. + */ struct HitInfo { - int hitid = -9; - int cryo = -9; - int tpc = -9; - int plane = -9; - int wire = -9; - int chan = -9; - float amp = -9; - float rms = -9; - int trkid = -9; - int shwrid = -9; - int clustid = -9; - int blipid = -9; - bool ismatch = false; - float integralADC = -999; // [ADCs] from integral - float sigmaintegral = -999; - float sumADC = -999; // [ADCs] from sum - float charge = -999; // [e-] - float peakTime = -999999; - float driftTime = -999999; // [tick] - float gof = -9; - int g4trkid = -9; - int g4pdg = -999; - int g4charge = -999; // [e-] - float g4frac = -99; - float g4energy = -999; // [MeV] - }; - - struct HitClust { - int ID = -9; - bool isValid = false; - int CenterChan = -999; - int CenterWire = -999; - bool isTruthMatched = false; - bool isMerged = false; - bool isMatched = false; - int DeadWireSep = 99; - int Cryostat = -9; - int TPC = -9; - int Plane = -9; - int NHits = -9; - int NWires = -9; - float ADCs = -999; - float Amplitude = -999; - float Charge = -999; - float SigmaCharge = -999; - float TimeTick = -999; - float Time = -999; - float StartHitTime = -999; - float EndHitTime = -999; - float StartTime = -999; - float EndTime = -999; - float Timespan = -999; - float RMS = -999; - int StartWire = -999; - int EndWire = -999; - int NPulseTrainHits = -9; - float GoodnessOfFit = -999; - int BlipID = -9; - int EdepID = -9; - si_t HitIDs; - si_t Wires; - si_t Chans; - si_t G4IDs; - - std::map IntersectLocations; + int hitid = -9; ///< Index of hit in recob::hit vector + int cryo = -9; ///< cryostat containing hit. From geo::WireID + int tpc = -9; ///< TPC containing hit. From geo::WireID + int plane = -9; ///< Plane containing hit. From geo::WireID + int wire = -9; ///< Wire containing hit. From geo::WireID + int chan = -9; ///< Channel ID containing hit. From recob::Hit->Channel() + float amp = -9; ///< Amplitude of hit [ADC] + float rms = -9; ///< RMS of hit shape [ticks] + int trkid = -9; ///< track ID from hit-track association object + int shwrid = -9; ///< Depreciated + int clustid = -9; ///< If this hit gets gathered into a cluster, what is the index of that cluster. All hits should be in a track or a cluster + int blipid = -9; ///< If this hit gets gathered into a blip, what is the index of that blip. + bool ismatch = false; ///< Depreciated + float integralADC = -999; ///< Integral area from the recob::hit [ADCs-ticks] + float sigmaintegral = -999; ///< Uncertainty on the integralADC value [ADC-ticks] + float sumADC = -999; ///< [ADCs] from ROISummedADC method + float charge = -999; ///< Total charge as estimated from hit integral [e-] + /*! + Charge is reconstructed from calo::CalorimetryAlg ( "sbnd_calorimetryalgmc" )-> ElectronsFromADCArea function + Configuration is in sbndcode/sbndcode/LArSoftConfigurations/calorimetry_sbnd.fcl + */ + float peakTime = -999999; ///< Hit peak time + a fcl-set per-plane offset [tick] + float driftTime = -999999; ///< Drift-time accounting for the interplane separation corrections [tick] + float gof = -9; ///< GoodnessOfFit from recob::hit + int g4trkid = -9; ///< simb::MCParticle g4 track ID associated with this hit + int g4pdg = -999; ///< simb::MCParticle g4pdg + int g4charge = -999; ///< anab::BackTrackerHitMatchingData numelectrons [e-] + float g4frac = -99; ///< anab::BackTrackerHitMatchingData ideNFraction, or fraction of energy in hit from this particle + float g4energy = -999; ///< anab::BackTrackerHitMatchingData energy [MeV] }; - struct Blip { - - int ID = -9; // Blip ID / index - bool isValid = false; // Blip passes basic checks - int Cryostat = -9; // Cryostat - int TPC = -9; // TPC - int NPlanes = -9; // Num. matched planes - int MaxWireSpan = -9; // Maximum span of wires on any plane cluster - float TimeTick = -999; // Readout time [ticks] - float Time = -999; // Drift time [us] - float Charge = -9; // Charge on calorimetry plane - float Energy = -999; // Energy (const dE/dx, fcl-configurable) - float EnergyESTAR = -999; // Energy (ESTAR method from ArgoNeuT) - float EnergyPSTAR = -999; // Energy (PSTAR method similar with ESTAR method from ArgoNeuT) - float ProxTrkDist = -9; // Distance to cloest track - int ProxTrkID = -9; // ID of closest track - bool inCylinder = false; // Is it in a cone/cylinder region? - - TVector3 Position; // 3D position TVector3 - float SigmaYZ = -9.; // Uncertainty in YZ intersect [cm] - float dX = -9; // Equivalent length along drift direction [cm] - float dYZ = -9; // Approximate length scale in YZ space [cm] +/** True energy depositions + * std::vector > makes a particle list + * std::vector makes a IDE list + * Together those make a blip::ParticleInfo object used to fill in this TrueBlip information + * For a reconstructed MC blip we will record the true blip info associated with it + * That blip reconstruction applies cuts to overall blip size/spread + * A single TrueBlip will be constructed for energy depositions within TrueBlipMergeDist (fcl set 0.3 cm by default) + */ + struct TrueBlip { + int ID = -9; ///< Index of this trueBlip object + int Cryostat = -9; ///< Cryostat ID the blip reconstructed to + int TPC = -9; ///< TPC ID the blip reconstructed to + float Time = -9; ///< Charge weighted peak time of TrueBlip energy depositions [tick] + int TimeTick = -9; ///< Depreciated + float DriftTime = -9; ///< Charge weighted drift time of TrueBlip energy depositions [tick] + float Energy = 0; ///< Total energy dep [MeV] + int DepElectrons = 0; ///< Total deposited electrons [e-] + int NumElectrons = 0; ///< electrons reaching wires [e-] + int LeadG4ID = -9; ///< G4 track ID depositing the most charge in this deposition + int LeadG4Index = -9; ///< G4 track Index depositing the most charge in this deposition + int LeadG4PDG = -9; ///< G4 PDG associated with the track depositing the most charge + float LeadCharge = -9; ///< Largest charge deposition associated with this True Blip + geo::Point_t Position; ///< Charge weighted true-XYZ position [cm] //Used to be TVector3 + typedef std::map G4ChargeMap; ///< Map from G4 particle track ID to deposited charge in this energy deposition + typedef std::map G4PDGMap; ///< Map from G4 particle track ID to G4 PDG + }; - // Plane/cluster-specific information - blip::HitClust clusters[kNplanes]; - - // Truth-matched energy deposition - blip::TrueBlip truth; - - // Prototype getter functions - double X() { return Position.X(); } - double Y() { return Position.Y(); } - double Z() { return Position.Z(); } - + /** Holder for MCTruth article data + * Filled from std::vector > (makes a particle list) + * and std::vector (makes a IDE list) + */ + struct ParticleInfo { + simb::MCParticle particle; + int trackId = -9; ///< ID from GEANT4 track. + int index = -9; ///< + int isPrimary = -9; ///< + int numTrajPts = -9; ///< + double depEnergy = -9; ///< + int depElectrons = -9; ///< + double numElectrons = -9; ///< + double mass = -9; ///< simb::MCParticle mass [MeV] + double E = -9; ///< Starting simb::MCParticle energy [MeV] + double endE = -9; ///< simb::MCParticle final energy [MeV] + double KE = -9; ///< E - Mass [MeV] + double endKE = -9; ///< EndE - Mass [MeV] + double P = -9; ///< Starting Magnitude of momentum vector [MeV/c] + double Px = -9; ///< X-component of the momentum vector [MeV/c] + double Py = -9; ///< Y-component of the momentum vector [MeV/c] + double Pz = -9; ///< Z-component of the momentum vector [MeV/c] + double pathLength = -9; ///< Sum of distances between interaction vertecies within the active volume [cm] + double time = -9; ///< Time at birth [us] + double endtime = -9; ///< Time of final interaction [us] + geo::Point_t startPoint; ///< Starting Position of particle [cm] //Used to be TVector3 + geo::Point_t endPoint; ///< End position of particle [cm] //Used to be TVector3 + geo::Point_t position; ///< Central position of the trajectory [cm] //Used to be TVector3 }; - } #endif diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index 60536524..941ed5b5 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -4,6 +4,7 @@ cet_make_library( LIBRARY_NAME sbndobj_BlipDataTypes LIBRARIES lardataobj::RecoBase nusimdata::SimulationBase + larcorealg::Geometry ) art_dictionary(DICTIONARY_LIBRARIES sbndobj_BlipDataTypes) install_headers() From 8b3cb45b2967d4406cee7a00a152e575bdd23134 Mon Sep 17 00:00:00 2001 From: root Date: Sat, 29 Nov 2025 11:31:29 -0600 Subject: [PATCH 15/41] changed to default cmake library --- sbnobj/SBND/Blip/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index 941ed5b5..1fe622c2 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -1,4 +1,4 @@ -cet_make_library( LIBRARY_NAME sbndobj_BlipDataTypes +cet_make_library( SOURCE BlipDataTypes.cc LIBRARIES @@ -6,7 +6,7 @@ cet_make_library( LIBRARY_NAME sbndobj_BlipDataTypes nusimdata::SimulationBase larcorealg::Geometry ) -art_dictionary(DICTIONARY_LIBRARIES sbndobj_BlipDataTypes) +art_dictionary(DICTIONARY_LIBRARIES sbnobj::SBND_Blip) install_headers() install_source() install_fhicl() From 47f8630a231e266984e89c709a5d99b39b5bcd68 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 11:47:47 -0600 Subject: [PATCH 16/41] Updating order in classes.h and adding class version index --- sbnobj/SBND/Blip/classes.h | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index 7508db02..0fbd2909 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -24,13 +24,17 @@ // Only include objects that we would like to be able to put into the event. // Do not include the objects they contain internally. // - -template class art::Assns; -template class art::Wrapper >; -template class std::vector; -template class art::Wrapper >; -template class std::map; -template class art::Assns; -template class art::Wrapper >; -template class art::Assns; -template class art::Wrapper >; +; +; +; +; +; +; +; +; +; +; +; +; +; +; \ No newline at end of file From d4991e582a2deffca59a767b076237952ffcc2c2 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 11:52:18 -0600 Subject: [PATCH 17/41] These have to be in a specific order --- sbnobj/SBND/Blip/BlipDataTypes.h | 199 ++++++++++++++++--------------- 1 file changed, 100 insertions(+), 99 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 94d1e546..1c0944d9 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -10,105 +10,6 @@ namespace blip { const int kNplanes = 3; ///< Wire planes in a given tpc in this detector const int kNTPCs = 2; ///< TPC in this detector - -/** Blips are small, plane-matched, energy deposition in liquid argon with sizes similar to wire separation. - * Blips are composed of hit-clusters, which are (time and space) adjacent hits on a single wire plane. - * A single plane, usually the collection, is given a privledged position as the calorimetry plane. - * Every hit-cluster on the calorimetry plane is checked for matches on the other two planes, and when multiple possible pairs are found - * the highest scoring one is selected. - * Score components include cluster charge, cluster time duration, cluster peak time, and wire intersections. - */ - struct Blip { - int ID = -9; ///< Internal index for blip labelling - bool isValid = false; ///< Blip passes basic checks in reco. All blips saved to artRoot file should be valid. - int Cryostat = -9; ///< Which cryostat the blip was reconstructed to. For SBND this should always be 0. - int TPC = -9; ///< Which tpc the blip was reconstructed to. For SBND this may be 0 or 1. - int NPlanes = -9; ///< Number of planes matched to build the blip. Blips must be plane matched so this should be 2+ planes. - int MaxWireSpan = -9; ///< Maximum span of wires on any plane hit-cluster - float TimeTick = -999; ///< Average time of hitclusters making up blips. [tpc tick] - /*! - Hit cluster time is the charge-weighted average of the hit-peak times for hits in the cluster - */ - float Time = -999; ///< Average time of hitclusters making up blips. [us] - /*! - Hit cluster time is the charge-weighted average of the hit-peak times for hits in the cluster - */ - float Charge = -9; ///< Charge on calorimetry plane [e-] - float Energy = -999; ///< Reconstructed energy in the calorimetry plane (const dE/dx, fcl-configurable) [MeV] - float EnergyESTAR = -999; ///< Reconstructed energy in the calorimetry plane (ESTAR method from ArgoNeuT) [MeV] - float EnergyPSTAR = -999; ///< Reconstructed energy in the calorimetry plane (PSTAR method similar with ESTAR method from ArgoNeuT) [MeV] - float ProxTrkDist = -9; ///< 3-D distance to closest track, assuming the blip was concident with event trigger [cm] - /*! - for properly flash matched out-of-time tracks this distance will be wrong! The blips have no such flash matching ability as of yet - */ - int ProxTrkID = -9; ///< index of the of closest track, assuming the blip was concident with event trigger [cm] - bool inCylinder = false; ///< Bool for whether this blip falls in a 45 degree cone relative to any track closer than fcl-set CylinderRadius (15 cm) - /*! - please note the blip X position is unreliable, so these distance and 3-d position derived variables may be incorrect - */ - geo::Point_t Position; ///< 3D position TVector3 //Used to be TVector3 - float SigmaYZ = -9.; ///< Uncertainty in YZ intersect [cm] - float dX = -9; ///< Equivalent length along drift direction [cm] - float dYZ = -9; ///< Approximate length scale in YZ space [cm] - blip::HitClust clusters[kNplanes]; ///< Plane/cluster-specific information. Just contains hit clusters making up this blip! - blip::TrueBlip truth; ///< Truth-matched energy deposition. Just contains the relevant MC truth info to this blip! - double X() { return Position.X(); } - double Y() { return Position.Y(); } - double Z() { return Position.Z(); } - - }; - - /** Hit clusters are collections of adjacent hits on a single wire plane. - * They have a fcl-set maximum wire-span (fcl set. Default 10), as well as maximum tick-width (fcl set. Default 50 tick). - * Hit clusters have timestamp and associated wire IDs. - * Within a hit cluster certain statistical summaries of the collection are saved including: - * Total charge, total charge uncertianty, peak hit amplitude, and charge weighted RMS hit spread - */ - struct HitClust { - int ID = -9; ///< Per-plane index for the hit clusters. In SBND we save every collection plane hitcluster but not the induction - bool isValid = false; ///< Bool check that every hit is in the same cryostat, tpc, plane. Should always be true for saved items - int CenterChan = -999; ///< Channel ID of the wire in the geometric center of the hit cluster - int CenterWire = -999; ///< Wire ID of the wire in the geometric center of the hit cluster - bool isTruthMatched = false; ///< is there a trueBlip with the same leadG4ID as one of the G4IDs making up this cluster - bool isMerged = false; ///< Depreciated internal flag for tracking blips before/after merge. - bool isMatched = false; ///< Is this hit cluster plane-matched into a full 3d blip - int DeadWireSep = 99; ///< Separation between the extreme ends of the hitcluster and the nearest dead wire. - /*! - DeadWireSep can be between 0 and 5 and valid. Larger separations are filled in as 99. - */ - int Cryostat = -9; ///< cryostat for this hit cluster - int TPC = -9; ///< TPC for this hit cluster - int Plane = -9; ///< Plane index for this hit cluster - int NHits = -9; ///< Number of hits making up this hit cluster - int NWires = -9; ///< Wire span of the hit cluster - float ADCs = -999; ///< ADC integral sum of hits making up the cluster [ADC-tick] - float Amplitude = -999; ///< Max amplitude of hits making up the hit cluster [ADC] - float Charge = -999; ///< Total charge of hits making up the hit cluster [e-] - /*! - Charge is reconstructed from calo::CalorimetryAlg ( "sbnd_calorimetryalgmc" )-> ElectronsFromADCArea function - Configuration is in sbndcode/sbndcode/LArSoftConfigurations/calorimetry_sbnd.fcl - */ - float SigmaCharge = -999; ///< charge-weighted charge uncertainties for this hit-cluster [e-] - float TimeTick = -999; ///< charge-weighted average hit-peak-times for this hit-cluster [tick] - float Time = -999; ///< charge-weighted average hit-peak-times for this hit-cluster [us] - float StartHitTime = -999; ///< Depreciated - float EndHitTime = -999; ///< Depreciated - float StartTime = -999; ///< Minimum -1 sigma time of a hit in this cluster [us] - float EndTime = -999; ///< Max +1 sigma time of a hit in this cluster [us] - float Timespan = -999; ///< Hit cluster EndTime - StartTime [us] - float RMS = -999; ///< Quadrature estimate of charge waveform timespread accounting for varied hit-drift times and internal hit-RMS [us] - int StartWire = -999; ///< Lowest wireID involved with the hitcluster - int EndWire = -999; ///< Highest wireID involved with the hit cluster - int NPulseTrainHits = -9; ///< Number of hits with a GoodnessOfFit<0 involved in this hit cluster - float GoodnessOfFit = -999; ///< Charge weighted hit-GoodnessofFit param - int BlipID = -9; ///< If this hit cluster ended up in a blip, what is its ID - int EdepID = -9; ///< If this hit cluster is MC-matched what is the trueBlip ID - std::set HitIDs; ///< Index of the recob::hit objects making up this cluster - std::set Wires; ///< Set of geo::wireIDs contributing hits to this cluster - std::set Chans; ///< Set of raw::ChannelID_t contributing hits to this cluster - std::set G4IDs; ///< simb::MCParticle track ID contributing hits to this cluster - std::map IntersectLocations; ///< Internal reconstruction variable for recording where hit-clusters on different planes would overlap //used to be TVector3 - }; /** Extra struct used in blip reco for grabbing key recob::hit information * Most attributes are grabbed directly from the recob::hit, but some, like charge, involve extra processing. @@ -201,5 +102,105 @@ namespace blip { geo::Point_t endPoint; ///< End position of particle [cm] //Used to be TVector3 geo::Point_t position; ///< Central position of the trajectory [cm] //Used to be TVector3 }; + + /** Hit clusters are collections of adjacent hits on a single wire plane. + * They have a fcl-set maximum wire-span (fcl set. Default 10), as well as maximum tick-width (fcl set. Default 50 tick). + * Hit clusters have timestamp and associated wire IDs. + * Within a hit cluster certain statistical summaries of the collection are saved including: + * Total charge, total charge uncertianty, peak hit amplitude, and charge weighted RMS hit spread + */ + struct HitClust { + int ID = -9; ///< Per-plane index for the hit clusters. In SBND we save every collection plane hitcluster but not the induction + bool isValid = false; ///< Bool check that every hit is in the same cryostat, tpc, plane. Should always be true for saved items + int CenterChan = -999; ///< Channel ID of the wire in the geometric center of the hit cluster + int CenterWire = -999; ///< Wire ID of the wire in the geometric center of the hit cluster + bool isTruthMatched = false; ///< is there a trueBlip with the same leadG4ID as one of the G4IDs making up this cluster + bool isMerged = false; ///< Depreciated internal flag for tracking blips before/after merge. + bool isMatched = false; ///< Is this hit cluster plane-matched into a full 3d blip + int DeadWireSep = 99; ///< Separation between the extreme ends of the hitcluster and the nearest dead wire. + /*! + DeadWireSep can be between 0 and 5 and valid. Larger separations are filled in as 99. + */ + int Cryostat = -9; ///< cryostat for this hit cluster + int TPC = -9; ///< TPC for this hit cluster + int Plane = -9; ///< Plane index for this hit cluster + int NHits = -9; ///< Number of hits making up this hit cluster + int NWires = -9; ///< Wire span of the hit cluster + float ADCs = -999; ///< ADC integral sum of hits making up the cluster [ADC-tick] + float Amplitude = -999; ///< Max amplitude of hits making up the hit cluster [ADC] + float Charge = -999; ///< Total charge of hits making up the hit cluster [e-] + /*! + Charge is reconstructed from calo::CalorimetryAlg ( "sbnd_calorimetryalgmc" )-> ElectronsFromADCArea function + Configuration is in sbndcode/sbndcode/LArSoftConfigurations/calorimetry_sbnd.fcl + */ + float SigmaCharge = -999; ///< charge-weighted charge uncertainties for this hit-cluster [e-] + float TimeTick = -999; ///< charge-weighted average hit-peak-times for this hit-cluster [tick] + float Time = -999; ///< charge-weighted average hit-peak-times for this hit-cluster [us] + float StartHitTime = -999; ///< Depreciated + float EndHitTime = -999; ///< Depreciated + float StartTime = -999; ///< Minimum -1 sigma time of a hit in this cluster [us] + float EndTime = -999; ///< Max +1 sigma time of a hit in this cluster [us] + float Timespan = -999; ///< Hit cluster EndTime - StartTime [us] + float RMS = -999; ///< Quadrature estimate of charge waveform timespread accounting for varied hit-drift times and internal hit-RMS [us] + int StartWire = -999; ///< Lowest wireID involved with the hitcluster + int EndWire = -999; ///< Highest wireID involved with the hit cluster + int NPulseTrainHits = -9; ///< Number of hits with a GoodnessOfFit<0 involved in this hit cluster + float GoodnessOfFit = -999; ///< Charge weighted hit-GoodnessofFit param + int BlipID = -9; ///< If this hit cluster ended up in a blip, what is its ID + int EdepID = -9; ///< If this hit cluster is MC-matched what is the trueBlip ID + std::set HitIDs; ///< Index of the recob::hit objects making up this cluster + std::set Wires; ///< Set of geo::wireIDs contributing hits to this cluster + std::set Chans; ///< Set of raw::ChannelID_t contributing hits to this cluster + std::set G4IDs; ///< simb::MCParticle track ID contributing hits to this cluster + std::map IntersectLocations; ///< Internal reconstruction variable for recording where hit-clusters on different planes would overlap //used to be TVector3 + }; + + + /** Blips are small, plane-matched, energy deposition in liquid argon with sizes similar to wire separation. + * Blips are composed of hit-clusters, which are (time and space) adjacent hits on a single wire plane. + * A single plane, usually the collection, is given a privledged position as the calorimetry plane. + * Every hit-cluster on the calorimetry plane is checked for matches on the other two planes, and when multiple possible pairs are found + * the highest scoring one is selected. + * Score components include cluster charge, cluster time duration, cluster peak time, and wire intersections. + */ + struct Blip { + int ID = -9; ///< Internal index for blip labelling + bool isValid = false; ///< Blip passes basic checks in reco. All blips saved to artRoot file should be valid. + int Cryostat = -9; ///< Which cryostat the blip was reconstructed to. For SBND this should always be 0. + int TPC = -9; ///< Which tpc the blip was reconstructed to. For SBND this may be 0 or 1. + int NPlanes = -9; ///< Number of planes matched to build the blip. Blips must be plane matched so this should be 2+ planes. + int MaxWireSpan = -9; ///< Maximum span of wires on any plane hit-cluster + float TimeTick = -999; ///< Average time of hitclusters making up blips. [tpc tick] + /*! + Hit cluster time is the charge-weighted average of the hit-peak times for hits in the cluster + */ + float Time = -999; ///< Average time of hitclusters making up blips. [us] + /*! + Hit cluster time is the charge-weighted average of the hit-peak times for hits in the cluster + */ + float Charge = -9; ///< Charge on calorimetry plane [e-] + float Energy = -999; ///< Reconstructed energy in the calorimetry plane (const dE/dx, fcl-configurable) [MeV] + float EnergyESTAR = -999; ///< Reconstructed energy in the calorimetry plane (ESTAR method from ArgoNeuT) [MeV] + float EnergyPSTAR = -999; ///< Reconstructed energy in the calorimetry plane (PSTAR method similar with ESTAR method from ArgoNeuT) [MeV] + float ProxTrkDist = -9; ///< 3-D distance to closest track, assuming the blip was concident with event trigger [cm] + /*! + for properly flash matched out-of-time tracks this distance will be wrong! The blips have no such flash matching ability as of yet + */ + int ProxTrkID = -9; ///< index of the of closest track, assuming the blip was concident with event trigger [cm] + bool inCylinder = false; ///< Bool for whether this blip falls in a 45 degree cone relative to any track closer than fcl-set CylinderRadius (15 cm) + /*! + please note the blip X position is unreliable, so these distance and 3-d position derived variables may be incorrect + */ + geo::Point_t Position; ///< 3D position TVector3 //Used to be TVector3 + float SigmaYZ = -9.; ///< Uncertainty in YZ intersect [cm] + float dX = -9; ///< Equivalent length along drift direction [cm] + float dYZ = -9; ///< Approximate length scale in YZ space [cm] + blip::HitClust clusters[kNplanes]; ///< Plane/cluster-specific information. Just contains hit clusters making up this blip! + blip::TrueBlip truth; ///< Truth-matched energy deposition. Just contains the relevant MC truth info to this blip! + double X() { return Position.X(); } + double Y() { return Position.Y(); } + double Z() { return Position.Z(); } + + }; } #endif From 91eba8c8a997a435461115a61d1d88d1a65db064 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 12:00:37 -0600 Subject: [PATCH 18/41] Getting rid of ; --- sbnobj/SBND/Blip/classes.h | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index 0fbd2909..ece98e96 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -24,17 +24,17 @@ // Only include objects that we would like to be able to put into the event. // Do not include the objects they contain internally. // -; -; -; -; -; -; -; -; -; -; -; -; -; -; \ No newline at end of file + + + + + + + + + + + + + + \ No newline at end of file From 83efc1025df2b04059bbf15b045f1e9f0d180c67 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 13:05:26 -0600 Subject: [PATCH 19/41] don't need this one --- sbnobj/SBND/Blip/classes.h | 1 - 1 file changed, 1 deletion(-) diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index ece98e96..93d2af9a 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -24,7 +24,6 @@ // Only include objects that we would like to be able to put into the event. // Do not include the objects they contain internally. // - From 3d4c322be2fdf61db1e7f38cbd090770b56a306c Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 13:11:28 -0600 Subject: [PATCH 20/41] Adding blank lines to be filled by compiler --- sbnobj/SBND/Blip/classes.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index 93d2af9a..98932c9e 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -25,12 +25,15 @@ // Do not include the objects they contain internally. // + + + From 4e06944016ddaf0bdd3c105baac08cbc3ab2d314 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 13:13:48 -0600 Subject: [PATCH 21/41] Adding blank lines to be filled by compiler --- sbnobj/SBND/Blip/classes.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index 98932c9e..fb906c76 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -25,15 +25,15 @@ // Do not include the objects they contain internally. // - + - + - + From 5ae8d148f027ade10441d7b5f9d1a2c6f495f096 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 13:16:27 -0600 Subject: [PATCH 22/41] I don't understand this version tracker --- sbnobj/SBND/Blip/classes.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index fb906c76..037fb596 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -24,15 +24,15 @@ // Only include objects that we would like to be able to put into the event. // Do not include the objects they contain internally. // - + - + - + From 1fd419a42029b66ade81a2db402bff8277e042d3 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 13:20:59 -0600 Subject: [PATCH 23/41] I was putting things in the wrong place --- sbnobj/SBND/Blip/classes.h | 26 ++++++++++---------------- sbnobj/SBND/Blip/classes_def.xml | 8 ++++---- 2 files changed, 14 insertions(+), 20 deletions(-) diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index 037fb596..9e96d94e 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -24,19 +24,13 @@ // Only include objects that we would like to be able to put into the event. // Do not include the objects they contain internally. // - - - - - - - - - - - - - - - - \ No newline at end of file + +template class art::Assns; +template class art::Wrapper >; +template class std::vector; +template class art::Wrapper >; +template class std::map; +template class art::Assns; +template class art::Wrapper >; +template class art::Assns; +template class art::Wrapper >; diff --git a/sbnobj/SBND/Blip/classes_def.xml b/sbnobj/SBND/Blip/classes_def.xml index b6cf251f..28e2f91e 100644 --- a/sbnobj/SBND/Blip/classes_def.xml +++ b/sbnobj/SBND/Blip/classes_def.xml @@ -1,12 +1,12 @@ - + - + - - + + From 913d94bee3cec03bb87f0a5730729092138268f1 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 13:23:31 -0600 Subject: [PATCH 24/41] fixed stray typedef --- sbnobj/SBND/Blip/BlipDataTypes.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 1c0944d9..7d0a3b9d 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -69,8 +69,8 @@ namespace blip { int LeadG4PDG = -9; ///< G4 PDG associated with the track depositing the most charge float LeadCharge = -9; ///< Largest charge deposition associated with this True Blip geo::Point_t Position; ///< Charge weighted true-XYZ position [cm] //Used to be TVector3 - typedef std::map G4ChargeMap; ///< Map from G4 particle track ID to deposited charge in this energy deposition - typedef std::map G4PDGMap; ///< Map from G4 particle track ID to G4 PDG + std::map G4ChargeMap; ///< Map from G4 particle track ID to deposited charge in this energy deposition + std::map G4PDGMap; ///< Map from G4 particle track ID to G4 PDG }; /** Holder for MCTruth article data From ec53ab8808436b2fefadb8dd2455c9f31549ee22 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 13:31:56 -0600 Subject: [PATCH 25/41] Added checksum --- sbnobj/SBND/Blip/classes_def.xml | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/sbnobj/SBND/Blip/classes_def.xml b/sbnobj/SBND/Blip/classes_def.xml index 28e2f91e..9a0bcec5 100644 --- a/sbnobj/SBND/Blip/classes_def.xml +++ b/sbnobj/SBND/Blip/classes_def.xml @@ -1,11 +1,17 @@ - + + + - + + + - + + + From 1595e6cae3954a40957d7eaf16624d394d0dba8e Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 14:07:20 -0600 Subject: [PATCH 26/41] Missed a few comments --- sbnobj/SBND/Blip/BlipDataTypes.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 7d0a3b9d..cc76ccfb 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -80,12 +80,12 @@ namespace blip { struct ParticleInfo { simb::MCParticle particle; int trackId = -9; ///< ID from GEANT4 track. - int index = -9; ///< - int isPrimary = -9; ///< - int numTrajPts = -9; ///< - double depEnergy = -9; ///< - int depElectrons = -9; ///< - double numElectrons = -9; ///< + int index = -9; ///< Index in vector of simb::MCParticle + int isPrimary = -9; ///< Bool value of simbMCParticle is Primary + int numTrajPts = -9; ///< part.NumberTrajectoryPoints() from mcparticle + double depEnergy = -9; ///< Total energy deposited across all interactions from this partcle [MeV] + int depElectrons = -9; ///< Total electrons deposited across all interactions from this partcle [e-] + double numElectrons = -9; ///< Depreciated double mass = -9; ///< simb::MCParticle mass [MeV] double E = -9; ///< Starting simb::MCParticle energy [MeV] double endE = -9; ///< simb::MCParticle final energy [MeV] From 6194c0337c1e1dc55aea812f3b89960e3e004997 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 14:13:51 -0600 Subject: [PATCH 27/41] Changed array initilization --- sbnobj/SBND/Blip/BlipDataTypes.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index cc76ccfb..2adb93e8 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -195,7 +195,7 @@ namespace blip { float SigmaYZ = -9.; ///< Uncertainty in YZ intersect [cm] float dX = -9; ///< Equivalent length along drift direction [cm] float dYZ = -9; ///< Approximate length scale in YZ space [cm] - blip::HitClust clusters[kNplanes]; ///< Plane/cluster-specific information. Just contains hit clusters making up this blip! + std::array clusters; ///< Plane/cluster-specific information. Just contains hit clusters making up this blip! blip::TrueBlip truth; ///< Truth-matched energy deposition. Just contains the relevant MC truth info to this blip! double X() { return Position.X(); } double Y() { return Position.Y(); } From 41e66f4b7678f6d18e77960743ef19abcde0cae8 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 14:16:08 -0600 Subject: [PATCH 28/41] Removing some classes.h --- sbnobj/SBND/Blip/classes.h | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index 9e96d94e..bf163604 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -25,12 +25,11 @@ // Do not include the objects they contain internally. // -template class art::Assns; -template class art::Wrapper >; +template class art::Assns; +template class art::Wrapper >; template class std::vector; template class art::Wrapper >; -template class std::map; -template class art::Assns; -template class art::Wrapper >; -template class art::Assns; -template class art::Wrapper >; +template class art::Assns; +template class art::Wrapper >; +template class art::Assns; +template class art::Wrapper >; From 3736dddc3e48a1ce1806a1a4b118baf942833f15 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 14:16:32 -0600 Subject: [PATCH 29/41] Updated Checksum --- sbnobj/SBND/Blip/classes_def.xml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sbnobj/SBND/Blip/classes_def.xml b/sbnobj/SBND/Blip/classes_def.xml index 9a0bcec5..d670aa9b 100644 --- a/sbnobj/SBND/Blip/classes_def.xml +++ b/sbnobj/SBND/Blip/classes_def.xml @@ -1,6 +1,7 @@ - + + From cc2445271eabeb499b341bd3929bb2ce4eb5b727 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 14:18:40 -0600 Subject: [PATCH 30/41] Changing class version --- sbnobj/SBND/Blip/classes_def.xml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/sbnobj/SBND/Blip/classes_def.xml b/sbnobj/SBND/Blip/classes_def.xml index d670aa9b..ba196601 100644 --- a/sbnobj/SBND/Blip/classes_def.xml +++ b/sbnobj/SBND/Blip/classes_def.xml @@ -1,8 +1,7 @@ - - - + + From dfd4d1d7abb6a1eb6ac36d53bf84af564a9b08b4 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 14:19:52 -0600 Subject: [PATCH 31/41] remove unsaved item --- sbnobj/SBND/Blip/classes_def.xml | 1 - 1 file changed, 1 deletion(-) diff --git a/sbnobj/SBND/Blip/classes_def.xml b/sbnobj/SBND/Blip/classes_def.xml index ba196601..a8ef45ae 100644 --- a/sbnobj/SBND/Blip/classes_def.xml +++ b/sbnobj/SBND/Blip/classes_def.xml @@ -12,7 +12,6 @@ - From 4c9288e84308a2f288e0c97583426dd4686484bb Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sat, 29 Nov 2025 14:22:16 -0600 Subject: [PATCH 32/41] removed unused headers --- sbnobj/SBND/Blip/classes.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/sbnobj/SBND/Blip/classes.h b/sbnobj/SBND/Blip/classes.h index bf163604..5c902457 100644 --- a/sbnobj/SBND/Blip/classes.h +++ b/sbnobj/SBND/Blip/classes.h @@ -12,11 +12,8 @@ // data-products // lardataobj -//#include "lardata/Utilities/AssociationUtil.h" #include "canvas/Persistency/Common/Assns.h" -#include "lardataobj/RecoBase/PFParticle.h" #include "lardataobj/RecoBase/Hit.h" -#include "nusimdata/SimulationBase/MCTruth.h" #include "sbnobj/SBND/Blip/BlipDataTypes.h" #include "lardataobj/RecoBase/SpacePoint.h" From d4608537dc166adb3bacfd5cea53376814c7c809 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Sun, 30 Nov 2025 16:32:36 -0600 Subject: [PATCH 33/41] Fixed a comment --- sbnobj/SBND/Blip/BlipDataTypes.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 2adb93e8..29d04a17 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -186,12 +186,12 @@ namespace blip { /*! for properly flash matched out-of-time tracks this distance will be wrong! The blips have no such flash matching ability as of yet */ - int ProxTrkID = -9; ///< index of the of closest track, assuming the blip was concident with event trigger [cm] + int ProxTrkID = -9; ///< index of the of closest track, assuming the blip was concident with event trigger bool inCylinder = false; ///< Bool for whether this blip falls in a 45 degree cone relative to any track closer than fcl-set CylinderRadius (15 cm) /*! please note the blip X position is unreliable, so these distance and 3-d position derived variables may be incorrect */ - geo::Point_t Position; ///< 3D position TVector3 //Used to be TVector3 + geo::Point_t Position; ///< 3D position vector. Reconstructed with wrong t0! [cm] float SigmaYZ = -9.; ///< Uncertainty in YZ intersect [cm] float dX = -9; ///< Equivalent length along drift direction [cm] float dYZ = -9; ///< Approximate length scale in YZ space [cm] From c16b5bed85be2aac8b4e284e36d6103802301536 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Mon, 5 Jan 2026 11:23:19 -0600 Subject: [PATCH 34/41] made pinfo position have a capital like all other blip positions --- sbnobj/SBND/Blip/BlipDataTypes.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 29d04a17..07d87198 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -100,7 +100,7 @@ namespace blip { double endtime = -9; ///< Time of final interaction [us] geo::Point_t startPoint; ///< Starting Position of particle [cm] //Used to be TVector3 geo::Point_t endPoint; ///< End position of particle [cm] //Used to be TVector3 - geo::Point_t position; ///< Central position of the trajectory [cm] //Used to be TVector3 + geo::Point_t Position; ///< Central position of the trajectory [cm] //Used to be TVector3 }; /** Hit clusters are collections of adjacent hits on a single wire plane. From f4872e11c21ac175348ec8bf4fec24a5e7e3fa65 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Mon, 5 Jan 2026 13:00:29 -0600 Subject: [PATCH 35/41] Trying to fix runtime error --- sbnobj/SBND/Blip/BlipDataTypes.h | 1 + 1 file changed, 1 insertion(+) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 07d87198..197758c9 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -3,6 +3,7 @@ #include "lardataobj/RecoBase/Hit.h" #include "nusimdata/SimulationBase/MCParticle.h" #include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h" +#include "larcore/Geometry/Geometry.h" #include #include #include From e718325b4a071467ea1c08eb316bfe03fed7b3bd Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Mon, 5 Jan 2026 13:03:03 -0600 Subject: [PATCH 36/41] Trying to fix runtime error --- sbnobj/SBND/Blip/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index 1fe622c2..8a46e53b 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -4,7 +4,7 @@ cet_make_library( LIBRARIES lardataobj::RecoBase nusimdata::SimulationBase - larcorealg::Geometry + larcore::Geometry ) art_dictionary(DICTIONARY_LIBRARIES sbnobj::SBND_Blip) install_headers() From 693044314741f2907c6d255cc3780fcdb0c54476 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Mon, 5 Jan 2026 13:05:20 -0600 Subject: [PATCH 37/41] Trying to fix runtime error --- sbnobj/SBND/Blip/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index 8a46e53b..1fe622c2 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -4,7 +4,7 @@ cet_make_library( LIBRARIES lardataobj::RecoBase nusimdata::SimulationBase - larcore::Geometry + larcorealg::Geometry ) art_dictionary(DICTIONARY_LIBRARIES sbnobj::SBND_Blip) install_headers() From 5b3feeed953ee37a410095a0029d9a1cd1b9fc65 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Mon, 5 Jan 2026 13:07:24 -0600 Subject: [PATCH 38/41] Trying to fix runtime error --- sbnobj/SBND/Blip/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index 1fe622c2..4b84af96 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -5,6 +5,7 @@ cet_make_library( lardataobj::RecoBase nusimdata::SimulationBase larcorealg::Geometry + larcore::Geometry_Geometry_service ) art_dictionary(DICTIONARY_LIBRARIES sbnobj::SBND_Blip) install_headers() From 55f7c5da17fecce0c6fd18e11bbac556fbab7854 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Mon, 5 Jan 2026 16:37:42 -0600 Subject: [PATCH 39/41] removing the thing I added --- sbnobj/SBND/Blip/BlipDataTypes.h | 1 - 1 file changed, 1 deletion(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 197758c9..08577fd2 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -2,7 +2,6 @@ #define BLIPDATATYPE #include "lardataobj/RecoBase/Hit.h" #include "nusimdata/SimulationBase/MCParticle.h" -#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h" #include "larcore/Geometry/Geometry.h" #include #include From f7b4b84cdde32ba135cf16bd98ccce71fc5f14e6 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Mon, 5 Jan 2026 16:59:14 -0600 Subject: [PATCH 40/41] fixing run time error --- sbnobj/SBND/Blip/classes_def.xml | 1 + 1 file changed, 1 insertion(+) diff --git a/sbnobj/SBND/Blip/classes_def.xml b/sbnobj/SBND/Blip/classes_def.xml index a8ef45ae..b3b4c46e 100644 --- a/sbnobj/SBND/Blip/classes_def.xml +++ b/sbnobj/SBND/Blip/classes_def.xml @@ -13,6 +13,7 @@ + From 9524c9e7cee6291c2a81fb1179071f550d0f0ac5 Mon Sep 17 00:00:00 2001 From: Jacob McLaughlin Date: Tue, 6 Jan 2026 18:05:22 -0600 Subject: [PATCH 41/41] The geometry wasn't needed. Geo vectors was --- sbnobj/SBND/Blip/BlipDataTypes.h | 2 +- sbnobj/SBND/Blip/CMakeLists.txt | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/sbnobj/SBND/Blip/BlipDataTypes.h b/sbnobj/SBND/Blip/BlipDataTypes.h index 08577fd2..07d87198 100644 --- a/sbnobj/SBND/Blip/BlipDataTypes.h +++ b/sbnobj/SBND/Blip/BlipDataTypes.h @@ -2,7 +2,7 @@ #define BLIPDATATYPE #include "lardataobj/RecoBase/Hit.h" #include "nusimdata/SimulationBase/MCParticle.h" -#include "larcore/Geometry/Geometry.h" +#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h" #include #include #include diff --git a/sbnobj/SBND/Blip/CMakeLists.txt b/sbnobj/SBND/Blip/CMakeLists.txt index 4b84af96..fd893907 100644 --- a/sbnobj/SBND/Blip/CMakeLists.txt +++ b/sbnobj/SBND/Blip/CMakeLists.txt @@ -4,8 +4,6 @@ cet_make_library( LIBRARIES lardataobj::RecoBase nusimdata::SimulationBase - larcorealg::Geometry - larcore::Geometry_Geometry_service ) art_dictionary(DICTIONARY_LIBRARIES sbnobj::SBND_Blip) install_headers()