Feature/adding blip to caf #871
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change | ||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
@@ -55,7 +55,9 @@ namespace BlipUtils { | |||||||||||||||||||||||||||||
| pinfo.pathLength = PathLength( part, pinfo.startPoint, pinfo.endPoint); | ||||||||||||||||||||||||||||||
|
|
||||||||||||||||||||||||||||||
| // Central position of trajectory | ||||||||||||||||||||||||||||||
| 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.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.
Outdated
There was a problem hiding this comment.
This is the hardest one, being weighted. I would use the accumulator directly:
| 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")
Outdated
There was a problem hiding this comment.
You can directly use .R():
| 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(). 🤷
There was a problem hiding this comment.
As above:
| 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_iat 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();
There was a problem hiding this comment.
At least on this one we get the weights slightly wrong doing the suggested splitting of additions to the accumulators, since we have to divide by the total charge.
Since we don't know which one blip_j will be its hard to apply the right contribution in the outer loop.
To avoid overadding blip_i after the position update i call mpalg.clear();
There was a problem hiding this comment.
Scratch out the last comment. I couldn't get the compiler happy with the mpalg.add calls
In file included from /cvmfs/larsoft.opensciencegrid.org/products/lardataobj/v10_03_02/include/lardataobj/RecoBase/TrackingTypes.h:5, from /cvmfs/larsoft.opensciencegrid.org/products/lardataobj/v10_03_02/include/lardataobj/RecoBase/Trajectory.h:18, from /cvmfs/larsoft.opensciencegrid.org/products/lardataobj/v10_03_02/include/lardataobj/RecoBase/TrackTrajectory.h:19, from /cvmfs/larsoft.opensciencegrid.org/products/lardataobj/v10_03_02/include/lardataobj/RecoBase/Track.h:13, from /exp/sbnd/app/users/jmclaugh/sbndcode_v10_14_00_01/srcs/sbndcode/sbndcode/BlipRecoSBND/Utils/BlipUtils.h:23, from /exp/sbnd/app/users/jmclaugh/sbndcode_v10_14_00_01/srcs/sbndcode/sbndcode/BlipRecoSBND/Utils/BlipUtils.cc:1: /cvmfs/larsoft.opensciencegrid.org/products/larcorealg/v10_00_03/include/larcorealg/Geometry/geo_vectors_utils.h: In instantiation of ‘void geo::vect::MiddlePointAccumulatorDim<N>::add(BeginIter, EndIter) [with BeginIter = ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*; EndIter = float; unsigned int N = 3]’: /exp/sbnd/app/users/jmclaugh/sbndcode_v10_14_00_01/srcs/sbndcode/sbndcode/BlipRecoSBND/Utils/BlipUtils.cc:159:16: required from here /cvmfs/larsoft.opensciencegrid.org/products/larcorealg/v10_00_03/include/larcorealg/Geometry/geo_vectors_utils.h:1420:22: error: no matching function for call to ‘for_each(ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*&, float&, geo::vect::MiddlePointAccumulatorDim<3>::add<ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*, float>(ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*, float)::<lambda(const auto:16&)>)’ 1420 | std::for_each(begin, end, [this](auto const& p) { this->add(p); }); | ~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ In file included from /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/functional:64, from /cvmfs/larsoft.opensciencegrid.org/products/canvas/v3_16_04/include/canvas/Persistency/Provenance/ProductID.h:9, from /cvmfs/larsoft.opensciencegrid.org/products/canvas/v3_16_04/include/canvas/Persistency/Common/RefCore.h:5, from /cvmfs/larsoft.opensciencegrid.org/products/canvas/v3_16_04/include/canvas/Persistency/Common/Ptr.h:44, from /exp/sbnd/app/users/jmclaugh/sbndcode_v10_14_00_01/srcs/sbndcode/sbndcode/BlipRecoSBND/Utils/BlipUtils.h:12: /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/bits/stl_algo.h:3781:5: note: candidate: ‘template<class _IIter, class _Funct> _Funct std::for_each(_IIter, _IIter, _Funct)’ 3781 | for_each(_InputIterator __first, _InputIterator __last, _Function __f) | ^~~~~~~~ /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/bits/stl_algo.h:3781:5: note: template argument deduction/substitution failed: /cvmfs/larsoft.opensciencegrid.org/products/larcorealg/v10_00_03/include/larcorealg/Geometry/geo_vectors_utils.h:1420:22: note: deduced conflicting types for parameter ‘_IIter’ (‘ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*’ and ‘float’) 1420 | std::for_each(begin, end, [this](auto const& p) { this->add(p); }); | ~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ In file included from /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/algorithm:73, from /cvmfs/larsoft.opensciencegrid.org/products/cetlib/v3_18_02/slf7.x86_64.e26.prof/include/cetlib/container_algorithms.h:12, from /cvmfs/larsoft.opensciencegrid.org/products/canvas/v3_16_04/include/canvas/Persistency/Common/detail/aggregate.h:5, from /cvmfs/larsoft.opensciencegrid.org/products/canvas/v3_16_04/include/canvas/Persistency/Common/Wrapper.h:18, from /cvmfs/larsoft.opensciencegrid.org/products/canvas/v3_16_04/include/canvas/Persistency/Common/Ptr.h:45: /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/pstl/glue_algorithm_defs.h:42:1: note: candidate: ‘template<class _ExecutionPolicy, class _ForwardIterator, class _Function> __pstl::__internal::__enable_if_execution_policy<_ExecutionPolicy, void> std::for_each(_ExecutionPolicy&&, _ForwardIterator, _ForwardIterator, _Function)’ 42 | for_each(_ExecutionPolicy&& __exec, _ForwardIterator __first, _ForwardIterator __last, _Function __f); | ^~~~~~~~ /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/pstl/glue_algorithm_defs.h:42:1: note: template argument deduction/substitution failed: /cvmfs/larsoft.opensciencegrid.org/products/larcorealg/v10_00_03/include/larcorealg/Geometry/geo_vectors_utils.h:1420:22: note: deduced conflicting types for parameter ‘_ForwardIterator’ (‘float’ and ‘geo::vect::MiddlePointAccumulatorDim<3>::add<ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*, float>(ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*, float)::<lambda(const auto:16&)>’) 1420 | std::for_each(begin, end, [this](auto const& p) { this->add(p); }); | ~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ /cvmfs/larsoft.opensciencegrid.org/products/larcorealg/v10_00_03/include/larcorealg/Geometry/geo_vectors_utils.h: In instantiation of ‘void geo::vect::MiddlePointAccumulatorDim<N>::add(BeginIter, EndIter) [with BeginIter = const ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*; EndIter = float; unsigned int N = 3]’: /exp/sbnd/app/users/jmclaugh/sbndcode_v10_14_00_01/srcs/sbndcode/sbndcode/BlipRecoSBND/Utils/BlipUtils.cc:212:20: required from here /cvmfs/larsoft.opensciencegrid.org/products/larcorealg/v10_00_03/include/larcorealg/Geometry/geo_vectors_utils.h:1420:22: error: no matching function for call to ‘for_each(const ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*&, float&, geo::vect::MiddlePointAccumulatorDim<3>::add<const ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*, float>(const ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*, float)::<lambda(const auto:16&)>)’ /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/bits/stl_algo.h:3781:5: note: candidate: ‘template<class _IIter, class _Funct> _Funct std::for_each(_IIter, _IIter, _Funct)’ 3781 | for_each(_InputIterator __first, _InputIterator __last, _Function __f) | ^~~~~~~~ /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/bits/stl_algo.h:3781:5: note: template argument deduction/substitution failed: /cvmfs/larsoft.opensciencegrid.org/products/larcorealg/v10_00_03/include/larcorealg/Geometry/geo_vectors_utils.h:1420:22: note: deduced conflicting types for parameter ‘_IIter’ (‘const ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*’ and ‘float’) 1420 | std::for_each(begin, end, [this](auto const& p) { this->add(p); }); | ~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/pstl/glue_algorithm_defs.h:42:1: note: candidate: ‘template<class _ExecutionPolicy, class _ForwardIterator, class _Function> __pstl::__internal::__enable_if_execution_policy<_ExecutionPolicy, void> std::for_each(_ExecutionPolicy&&, _ForwardIterator, _ForwardIterator, _Function)’ 42 | for_each(_ExecutionPolicy&& __exec, _ForwardIterator __first, _ForwardIterator __last, _Function __f); | ^~~~~~~~ /cvmfs/larsoft.opensciencegrid.org/products/gcc/v12_1_0/Linux64bit+3.10-2.17/include/c++/12.1.0/pstl/glue_algorithm_defs.h:42:1: note: template argument deduction/substitution failed: /cvmfs/larsoft.opensciencegrid.org/products/larcorealg/v10_00_03/include/larcorealg/Geometry/geo_vectors_utils.h:1420:22: note: deduced conflicting types for parameter ‘_ForwardIterator’ (‘float’ and ‘geo::vect::MiddlePointAccumulatorDim<3>::add<const ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*, float>(const ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double>, ROOT::Math::GlobalCoordinateSystemTag>*, float)::<lambda(const auto:16&)>’) 1420 | std::for_each(begin, end, [this](auto const& p) { this->add(p); });
So I reverted these to simple loops not taking advantage of new functionality of point_t classes
There was a problem hiding this comment.
Sloppy programming by that utility author (i.e. myself).
The compiler saw add(tracking::Point_t, float) and it chose add(BeginIter, EndIter) over add(Point const&, double). Probably position.add(blip_i.Position, (double) weight_i) would have worked.
There was a problem hiding this comment.
This should work:
| 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()).
There was a problem hiding this comment.
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:(
#include "larcorealg/Geometry/geo_vectors_utils_TVector.h"); or using the more genericconvertTo():(for this one:
#include "larcorealg/Geometry/geo_vectors_utils.h").