Skip to content

Feature/adding blip to caf #871

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Jjm321814 wants to merge 40 commits into
base: develop
Choose a base branch
from
Open

Feature/adding blip to caf #871

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
40 commits
Select commit Hold shift + click to select a range
3abd95c
adjustments to new paths
Nov 6, 2025
4d2df24
Trying to get compile to solve
Nov 6, 2025
619cc9d
Assorted changes to things I don't actually use so the compiler would…
Nov 10, 2025
8c71d0b
updated sbnobj
Nov 10, 2025
bb7fca3
Forgot to take out cmake lists
Nov 13, 2025
04973a2
CMAKEList updates
Nov 13, 2025
a8816a2
Updating includes
Nov 13, 2025
7a1e130
Tried to separate the removal of classes.h but it doesn't actually work
Nov 13, 2025
70ee378
compiler is so fussy
Nov 14, 2025
9a41019
compile changes
Nov 14, 2025
8ab2400
merge issues
Nov 14, 2025
b2c99ca
updated caf maker fcl
Nov 14, 2025
2736872
fcl changes to include blips
Nov 14, 2025
dbcec4b
MC caf updates
Nov 14, 2025
f63553b
Fixing changes I didn't intend to make
Nov 19, 2025
064df71
Merge branch 'develop' into feature/AddingBlipToCAF
Nov 19, 2025
1fbb828
Fixed indentation
Jjm321814 Nov 20, 2025
1031a78
Had a mislabel on the fcl param
Nov 26, 2025
90533bc
Merge branch 'feature/BlipBugFix' into feature/AddingBlipToCAF
Nov 26, 2025
382de1d
fixed library name
Jjm321814 Nov 29, 2025
1808514
TVector3 Replacement updates
Jjm321814 Nov 29, 2025
d137623
more TVector 3 conversion
Jjm321814 Nov 29, 2025
184ba70
Update Power calls
Jjm321814 Nov 29, 2025
ae53aa5
Not sure this should function
Jjm321814 Nov 29, 2025
785cf59
One more cmake update
Jjm321814 Nov 29, 2025
a23b8fe
How did I miss so many
Jjm321814 Nov 29, 2025
311ee02
How did I miss so many
Jjm321814 Nov 29, 2025
a36a2b5
Weird these cmake issues appeared now and not earlier
Jjm321814 Nov 29, 2025
e7951b5
Adding includes that are probably unneeded
Jjm321814 Nov 29, 2025
1db5fd5
Oh I just updated one spot
Jjm321814 Nov 29, 2025
813ff32
Oh I just updated one spot
Jjm321814 Nov 29, 2025
d075348
I hope I don't just replace all these tvector 3
Jjm321814 Nov 29, 2025
cbcb056
updated middle induction plane hit amplitude to match Maria's work. F…
Jjm321814 Dec 4, 2025
d7ca349
Added two reco2 fcl that let us move back to single hitfinding stream…
Jjm321814 Dec 4, 2025
794d240
Being extra careful
Jjm321814 Dec 4, 2025
f886d99
fixing relabel in the new fcl
Jjm321814 Dec 4, 2025
3a30ea5
addressing Henry comments
Jjm321814 Dec 10, 2025
4b0b4c5
Merge branch 'feature/BlipBugFix' into feature/AddingBlipToCAF
Jjm321814 Dec 10, 2025
298f171
Tab align issue
Jjm321814 Dec 10, 2025
9c94f32
Tab align issue
Jjm321814 Dec 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1072,7 +1072,7 @@ namespace blip {
TVector3 p2(b.X(), b.Y(), b.Z() );
// TO-DO: if this track starts or ends at a TPC boundary,
// we should extend p1 or p2 to outside the AV to avoid blind spots
TVector3 bp = newBlip.Position;
TVector3 bp(newBlip.Position.X(), newBlip.Position.Y(), newBlip.Position.Z());
Copy link
Member

@PetrilloAtWork PetrilloAtWork Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Conversion geo::Point_t → TVector3:

Suggested change
TVector3 bp(newBlip.Position.X(), newBlip.Position.Y(), newBlip.Position.Z());
TVector3 bp = geo::vect::toTVector3(newBlip.Position);

(#include "larcorealg/Geometry/geo_vectors_utils_TVector.h"); or using the more generic convertTo():

Suggested change
TVector3 bp(newBlip.Position.X(), newBlip.Position.Y(), newBlip.Position.Z());
TVector3 bp = geo::vect::convertTo<TVector3>(newBlip.Position);

(for this one: #include "larcorealg/Geometry/geo_vectors_utils.h").

float d = BlipUtils::DistToLine(p1,p2,bp);
if( d > 0 ) {
// update closest trkdist
Expand Down
4 changes: 3 additions & 1 deletion sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@
#include "larcorealg/Geometry/GeometryCore.h"
#include "larreco/Calorimetry/CalorimetryAlg.h"
#include "art/Framework/Principal/Event.h"

#include "larcore/Geometry/WireReadout.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"

// Microboone includes
//#include "ubevt/Database/TPCEnergyCalib/TPCEnergyCalibService.h"
Expand All @@ -54,6 +55,7 @@

// Blip-specific utils
#include "sbndcode/BlipRecoSBND/Utils/BlipUtils.h"
#include "sbnobj/SBND/Blip/BlipDataTypes.h"

// ROOT stuff
#include "TH1D.h"
Expand Down
3 changes: 2 additions & 1 deletion sbndcode/BlipRecoSBND/Alg/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ art_make_library(LIBRARY_NAME sbndcode_BlipRecoAlg
lardataobj::RawData
lardataobj::RecoBase
lardata::RecoObjects
larreco::Calorimetry
larpandora::LArPandoraInterface
nusimdata::SimulationBase
cetlib::cetlib
Expand All @@ -33,8 +34,8 @@ art_make_library(LIBRARY_NAME sbndcode_BlipRecoAlg
ROOT::Gdml
sbndcode_CRTUtils
sbnobj::Common_CRT
#sbndcode_CosmicIdUtils
sbndcode_BlipUtils
sbnobj::SBND_Blip
)

install_headers()
Expand Down
4 changes: 3 additions & 1 deletion sbndcode/BlipRecoSBND/BlipAna_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
#include "lardataobj/AnalysisBase/Calorimetry.h"
#include "larevt/SpaceChargeServices/SpaceChargeService.h"
#include "cetlib/search_path.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"

// SBND-specific includes
#include "sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h"
Expand Down Expand Up @@ -75,7 +76,8 @@ const int kMaxTrks = 1000;
const int kMaxBlips = 5000;
const int kMaxG4 = 30000;
const int kMaxEDeps = 10000;
const int kMaxTrkPts = 2000;
const int kMaxTrkPts = 2000;
const int kNplanes = blip::kNplanes;

class BlipAna;

Expand Down
1 change: 1 addition & 0 deletions sbndcode/BlipRecoSBND/BlipRecoProducer_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@
#include "lardataobj/RawData/RawDigit.h"
#include "lardataobj/RawData/raw.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"

// C++ includes
#include <vector>
Expand Down
6 changes: 3 additions & 3 deletions sbndcode/BlipRecoSBND/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ set( MODULE_LIBRARIES larcorealg::Geometry
larevt::Filters
lardataobj::RawData
lardataobj::RecoBase
larreco::Calorimetry
larreco::Calorimetry
larreco::RecoAlg
lardata::RecoObjects
larpandora::LArPandoraInterface
Expand All @@ -34,12 +34,12 @@ set( MODULE_LIBRARIES larcorealg::Geometry
sbndcode_RecoUtils
sbndcode_OpDetSim
sbndcode_BlipUtils
sbnobj::SBND_Blip
sbndcode_BlipRecoAlg
)

add_subdirectory(Utils)
add_subdirectory(Alg)

add_subdirectory(Utils)
cet_build_plugin(BlipAna art::Module SOURCE BlipAna_module.cc LIBRARIES ${MODULE_LIBRARIES} )
cet_build_plugin(BlipRecoProducer art::Module SOURCE BlipRecoProducer_module.cc LIBRARIES ${MODULE_LIBRARIES} )

Expand Down
28 changes: 19 additions & 9 deletions sbndcode/BlipRecoSBND/Utils/BlipUtils.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,9 @@ namespace BlipUtils {
pinfo.pathLength = PathLength( part, pinfo.startPoint, pinfo.endPoint);

// Central position of trajectory
pinfo.position = 0.5*(pinfo.startPoint+pinfo.endPoint);
pinfo.position.SetXYZ(0.5*(pinfo.startPoint.X()+pinfo.endPoint.X()),
0.5*(pinfo.startPoint.Y()+pinfo.endPoint.Y()),
0.5*(pinfo.startPoint.Z()+pinfo.endPoint.Z()) );

Comment on lines +58 to 61
Copy link
Member

@PetrilloAtWork PetrilloAtWork Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is an unweighted middle point, so:

Suggested change
pinfo.position.SetXYZ(0.5*(pinfo.startPoint.X()+pinfo.endPoint.X()),
0.5*(pinfo.startPoint.Y()+pinfo.endPoint.Y()),
0.5*(pinfo.startPoint.Z()+pinfo.endPoint.Z()) );
pinfo.position = geo::vect::middlePoint({ pinfo.startPoint, pinfo.endPoint });

(#include "larcorealg/Geometry/geo_vector_utils.h")
Note the braces around the list of points (you can have any number of them).
This function is a simpler interface to MiddlePointAccumulator, which I'll mention below.

// Energy/charge deposited by this particle, found using SimEnergyDeposits
pinfo.depEnergy = 0;
Expand Down Expand Up @@ -155,7 +157,9 @@ namespace BlipUtils {
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.Position.SetXYZ( w1*tblip.Position.X() + w2*pinfo.position.X(),
w1*tblip.Position.Y() + w2*pinfo.position.Y(),
w1*tblip.Position.Z() + w2*pinfo.position.Z());
Comment on lines +160 to +162
Copy link
Member

@PetrilloAtWork PetrilloAtWork Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the hardest one, being weighted. I would use the accumulator directly:

Suggested change
tblip.Position.SetXYZ( w1*tblip.Position.X() + w2*pinfo.position.X(),
w1*tblip.Position.Y() + w2*pinfo.position.Y(),
w1*tblip.Position.Z() + w2*pinfo.position.Z());
geo::vect::MiddlePointAccumulator mpalg;
mpalg.add(tblip.Position, w1);
mpalg.add(pinfo.Position, w2);
tblip.Position = mpalg.middlePoint();

(#include "larcorealg/Geometry/geo_vector_utils.h")

tblip.Time = w1*tblip.Time + w2*pinfo.time;
tblip.LeadCharge = pinfo.depElectrons;
// ... if the particle isn't a match, show's over
Expand Down Expand Up @@ -196,15 +200,17 @@ namespace BlipUtils {
// check that the times are similar (we don't want to merge
// together a blip that happened much later but in the same spot)
if( fabs(blip_i.Time - blip_j.Time) > 5 ) continue;
float d = (blip_i.Position-blip_j.Position).Mag();
float d = TMath::Sqrt((blip_i.Position-blip_j.Position).Mag2());
Copy link
Member

@PetrilloAtWork PetrilloAtWork Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can directly use .R():

Suggested change
float d = TMath::Sqrt((blip_i.Position-blip_j.Position).Mag2());
float d = (blip_i.Position-blip_j.Position).R();

Also note that for determining the closest point, you don't really need to take the square root of the square distance, because the minimum distance is also the minimum distance square, so you can look for the latter and save a tiny bit of computing time.
An awkward feature of the new vectors is that they have Mag2() and R(), but not Mag() and R2(). 🤷

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.Position.SetXYZ( w1*blip_i.Position.X() + w2*blip_j.Position.X(),
w1*blip_i.Position.Y() + w2*blip_j.Position.Y(),
w1*blip_i.Position.Z() + w2*blip_j.Position.Z());
Comment on lines +211 to +213
Copy link
Member

@PetrilloAtWork PetrilloAtWork Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As above:

Suggested change
blip_i.Position.SetXYZ( w1*blip_i.Position.X() + w2*blip_j.Position.X(),
w1*blip_i.Position.Y() + w2*blip_j.Position.Y(),
w1*blip_i.Position.Z() + w2*blip_j.Position.Z());
geo::vect::MiddlePointAccumulator mpalg;
mpalg.add(blip_i.Position, w1);
mpalg.add(blip_j.Position, w2);
blip_i.Position = mpalg.middlePoint();

However, here we are in a loop context, so I recommend that the accumulation is done first, and then the results saved. You can declare your accumulators (MiddlePointAccumulator, StatCollector) in the outer loop, dealing with blip_i (and feeding in there the blip i), and then fill them in the inner loop:

  • initialization in outer loop 
        for(size_t i=0; i<vtb.size(); i++){
      if( isGrouped.at(i) ) continue;
      else isGrouped.at(i) = true;
      auto& blip_i = vtb.at(i);
      geo::vect::MiddlePointAccumulator position;
      lar::util::StatCollector<double> time;
      float const weight_i = blip_i.DepElectrons;
      position.add(blip_i.Position, weight_i);
      time.add(blip_i.time, weight_i);
      // ...
      for(size_t j=i+1; j<vtb.size(); j++){
  • filling in here:
    Suggested change
    blip_i.Position.SetXYZ( w1*blip_i.Position.X() + w2*blip_j.Position.X(),
    w1*blip_i.Position.Y() + w2*blip_j.Position.Y(),
    w1*blip_i.Position.Z() + w2*blip_j.Position.Z());
    float const weight_j = blip_j.DepElectrons;
    position.add(blip_j.Position.X(), weight_j);
    time.add(blip_j.time, weight_j);
    // ...
  • putting back into blip_i at the end of the loop: 
          }//loop over blip_j
      blip_i.Position = position.middlePoint();
      blip_i.time = time.Average();
      // ...
      blip_i.ID = vtb_merged.size();

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;
Expand Down Expand Up @@ -419,12 +425,16 @@ namespace BlipUtils {
// 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];
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;
for(auto& v : wirex ) newblip.Position.SetXYZ( newblip.Position.X() + v.X() * fact,
newblip.Position.Y() + v.Y() * fact,
newblip.Position.Z() + v.Z() * fact);
Comment on lines +432 to +434
Copy link
Member

@PetrilloAtWork PetrilloAtWork Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should work:

Suggested change
for(auto& v : wirex ) newblip.Position.SetXYZ( newblip.Position.X() + v.X() * fact,
newblip.Position.Y() + v.Y() * fact,
newblip.Position.Z() + v.Z() * fact);
geo::vect::MiddlePointAccumulator position;
position.add(begin(wirex), end(wirex));
newblip.Position = position.middlePoint();

I think this will work, still, given the following, it would be better having wirex redefined as std::vector<geo::Point_t>, especially considering that it is used only for this computation, and much better to have two accumulators initialised before the loop (a StatCollector2D sc for y and z — x is 0 —, and optionally a MiddlePointAccumulator for the position — but that one can be extracted by the other two), collect the statistics adding the coordinates in the loop, and finally computing the σ (sqrt(sc.VarianceX() + sc.VarianceY())) and average position ({ 0, sc.AverageX(), sc.AverageY() }). Also, the definition of σ feels off to me: it's the average of distances, while the usual RMS is the square root of the average of the squared distances (as in my example with VarianceX()/Y()).

for(auto& v : wirex ) newblip.SigmaYZ += TMath::Sqrt( pow(v.X()-newblip.Position.X(), 2) +
pow(v.Y()-newblip.Position.Y(), 2) +
pow(v.Z()-newblip.Position.Z(), 2)) * 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;
Expand Down Expand Up @@ -613,7 +623,7 @@ namespace BlipUtils {

//=============================================================================
// Length of particle trajectory
double PathLength(const simb::MCParticle& part, TVector3& start, TVector3& end)
double PathLength(const simb::MCParticle& part, geo::Point_t& start, geo::Point_t& end)
{
int n = part.NumberTrajectoryPoints();
if( n <= 1 ) return 0.;
Expand All @@ -633,7 +643,7 @@ namespace BlipUtils {
return L;
}
double PathLength(const simb::MCParticle& part){
TVector3 a,b;
geo::Point_t a,b;
return PathLength(part,a,b);
}

Expand Down
7 changes: 3 additions & 4 deletions sbndcode/BlipRecoSBND/Utils/BlipUtils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,14 +27,13 @@
#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"
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"

// c++
#include <vector>
#include <map>

#include "sbndcode/BlipRecoSBND/Utils/DataTypes.h"
#include "sbnobj/SBND/Blip/BlipDataTypes.h"
#include "TH1D.h"


Expand Down Expand Up @@ -76,7 +75,7 @@ namespace BlipUtils {
//void HitTruth(art::Ptr<recob::Hit> const&, int&, float&, float&, float&);
//si_t HitTruthIds( art::Ptr<recob::Hit> const&);
//bool G4IdToMCTruth( int const, art::Ptr<simb::MCTruth>&);
double PathLength(const simb::MCParticle&, TVector3&, TVector3&);
double PathLength(const simb::MCParticle&, geo::Point_t&, geo::Point_t&);
double PathLength(const simb::MCParticle&);
bool IsAncestorOf(int, int, bool);
double DistToBoundary(const recob::Track::Point_t&);
Expand Down
3 changes: 1 addition & 2 deletions sbndcode/BlipRecoSBND/Utils/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@ cet_make_library(
LIBRARY_NAME
sbndcode_BlipUtils
LIBRARIES
PUBLIC
larcorealg::Geometry
larcore::Geometry_Geometry_service
larsim::Simulation
Expand Down Expand Up @@ -38,7 +37,7 @@ cet_make_library(
ROOT::Gdml
sbndcode_CRTUtils
sbnobj::Common_CRT
#sbndcode_CosmicIdUtils
sbnobj::SBND_Blip
)

art_dictionary(DICTIONARY_LIBRARIES sbndcode_BlipUtils)
Expand Down
174 changes: 0 additions & 174 deletions sbndcode/BlipRecoSBND/Utils/DataTypes.h
View file
Open in desktop

This file was deleted.

Loading