123 <<
"------------------------------------------------------------------------------";
139 std::unique_ptr<std::vector<recob::Vertex>> vcol(
new std::vector<recob::Vertex>);
140 std::unique_ptr<art::Assns<recob::Vertex, recob::Hit>> vhassn(
142 std::unique_ptr<art::Assns<recob::Vertex, recob::Track>> vtassn(
144 std::unique_ptr<art::Assns<recob::Vertex, recob::Shower>> vsassn(
147 std::vector<recob::Track>
const& trkIn = *trackListHandle;
149 mf::LogInfo(
"PrimaryVertexFinder") <<
"number of tracks in this event = " << trkIn.size();
152 std::vector<recob::SpacePoint> startpoints_vec;
154 std::vector<TVector3> startvec;
157 std::vector<TVector3> endvec;
160 std::vector<TVector3> dircosvec;
163 std::vector<std::pair<art::Ptr<recob::Track>,
double>> trackpair;
167 for (
unsigned int i = 0; i < trkIn.size(); ++i) {
169 std::tie(start, end) = trkIn[i].Extent();
170 startXYZ.SetXYZ(start.X(), start.Y(), start.Z());
171 endXYZ.SetXYZ(end.X(), end.Y(), end.Z());
173 double length = (endXYZ - startXYZ).Mag();
178 for (
size_t i = 0; i < trackpair.size(); ++i) {
179 mf::LogInfo(
"PrimaryVertexFinder") <<
"track id is = " << (trackpair[i].first)->ID()
180 <<
" track length = " << (trackpair[i].second);
183 std::sort(trackpair.rbegin(), trackpair.rend(),
sort_pred2);
185 mf::LogInfo(
"PrimaryVertexFinder") <<
"AFTER SORTING ";
186 for (
size_t i = 0; i < trackpair.size(); ++i) {
187 mf::LogInfo(
"PrimaryVertexFinder") <<
"track id is = " << (trackpair[i].first)->ID()
188 <<
" track length = " << (trackpair[i].second);
201 for (
size_t j = 0; j < trackpair.size(); ++j) {
206 std::vector<recob::SpacePoint const*>
const& spacepoints = TrackSpacePoints.at(track.
key());
208 startXYZ = trackpair[j].first->Vertex<TVector3>();
209 endXYZ = trackpair[j].first->End<TVector3>();
210 dircosXYZ = trackpair[j].first->VertexDirection<TVector3>();
212 startvec.push_back(startXYZ);
213 endvec.push_back(endXYZ);
214 dircosvec.push_back(dircosXYZ);
216 mf::LogInfo(
"PrimaryVertexFinder") <<
"PrimaryVertexFinder got " << spacepoints.size()
217 <<
" 3D spacepoint(s) from Track3Dreco.cxx";
220 startpoints_vec.emplace_back(spacepoints[0]->XYZ(),
221 spacepoints[0]->ErrXYZ(),
222 spacepoints[0]->Chisq(),
223 startpoints_vec.size());
227 for (
size_t i = 0; i < startvec.size(); ++i) {
228 mf::LogInfo(
"PrimaryVertexFinder") <<
"Tvector3 start point SORTED = ";
231 for (
size_t i = 0; i < dircosvec.size(); ++i) {
232 mf::LogInfo(
"PrimaryVertexFinder") <<
"Tvector3 dir cos SORTED = ";
233 dircosvec[i].Print();
236 std::vector<std::vector<int>> vertex_collection_int;
237 std::vector<std::vector<TVector3>> vertexcand_vec;
239 for (
unsigned int i = 0; i < trackpair.size(); ++i) {
240 for (
unsigned int j = i + 1; j < trackpair.size(); ++j) {
242 <<
"distance between " << i <<
" and " << j <<
" = " 244 double GAMMA =
gammavalue(startvec[i], startvec[j], dircosvec[i], dircosvec[j]);
245 double ALPHA =
alphavalue(GAMMA, startvec[i], startvec[j], dircosvec[i], dircosvec[j]);
247 MinDist(ALPHA, GAMMA, startvec[i], startvec[j], dircosvec[i], dircosvec[j]);
249 <<
"alpha = " << ALPHA <<
" gamma = " << GAMMA <<
" MINIMUM DISTANCE = " << MINDIST;
254 mf::LogInfo(
"PrimaryVertexFinder") <<
"POINTS ON THE TRACKS ARE:: ";
264 std::vector<int> newvertex_int;
265 std::vector<TVector3> vertexcand;
266 newvertex_int.push_back(i);
267 newvertex_int.push_back(j);
268 vertex_collection_int.push_back(newvertex_int);
270 vertexcand.push_back(TRACK1POINT);
271 vertexcand.push_back(TRACK2POINT);
272 vertexcand_vec.push_back(vertexcand);
278 vertex_collection_int[index].push_back(i);
279 vertexcand_vec[index].push_back(TRACK1POINT);
282 vertex_collection_int[index].push_back(j);
283 vertexcand_vec[index].push_back(TRACK2POINT);
291 for (
size_t i = 0; i < trackpair.size(); ++i) {
294 std::vector<int> temp;
295 std::vector<TVector3> temp1;
297 temp1.push_back(startvec[i]);
298 vertex_collection_int.push_back(temp);
299 vertexcand_vec.push_back(temp1);
305 std::vector<size_t> vTrackIndices, vShowerIndices;
316 for (
size_t i = 0; i < vertex_collection_int.size(); ++i) {
322 itr < vertex_collection_int[i].end();
325 <<
"vector elements at index " << i <<
" are " << *itr
326 <<
"\ntrack original ID = " << (trackpair[*itr].first)->ID();
328 vTrackIndices.push_back(trackpair[*itr].first.key());
330 mf::LogInfo(
"PrimaryVertexFinder") <<
"------------";
333 itr < vertexcand_vec[i].end();
339 elemsize = vertexcand_vec[i].size();
342 double avgx = x / elemsize;
343 double avgy = y / elemsize;
344 double avgz = z / elemsize;
346 Double_t vtxcoord[3];
352 vcol->push_back(the3Dvertex);
354 if (!vTrackIndices.empty()) {
357 evt, *vtassn, vcol->size() - 1, vTrackIndices.begin(), vTrackIndices.end());
358 for (
size_t tIndex : vTrackIndices) {
359 std::vector<art::Ptr<recob::Hit>>
const&
hits = TrackHits.at(tIndex);
362 vTrackIndices.clear();
365 if (!vShowerIndices.empty()) {
368 evt, *vsassn, vcol->size() - 1, vShowerIndices.begin(), vShowerIndices.end());
369 for (
size_t sIndex : vShowerIndices) {
370 std::vector<art::Ptr<recob::Hit>>
const& hits = ShowerHits.at(sIndex);
373 vShowerIndices.clear();
378 MF_LOG_VERBATIM(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
380 for (
size_t i = 0; i < vcol->size(); ++i)
383 evt.
put(std::move(vcol));
384 evt.
put(std::move(vtassn));
385 evt.
put(std::move(vhassn));
386 evt.
put(std::move(vsassn));
TVector3 PointOnExtendedTrack(double alphagamma, TVector3 startpoint, TVector3 dircos)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition of vertex object for LArSoft.
bool IsInNewVertex(int a, std::vector< int > newvertex)
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
key_type key() const noexcept
std::string fTrackModuleLabel
double MinDist(double alpha, double gamma, TVector3 startpoint1, TVector3 startpoint2, TVector3 dircos1, TVector3 dircos2)
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
double StartPointSeperation(recob::SpacePoint sp1, recob::SpacePoint sp2)
bool sort_pred2(const std::pair< art::Ptr< recob::Track >, double > &left, const std::pair< art::Ptr< recob::Track >, double > &right)
#define MF_LOG_VERBATIM(category)
tracking::Point_t Point_t
EventNumber_t event() const
double gammavalue(TVector3 startpoint1, TVector3 startpoint2, TVector3 dircos1, TVector3 dircos2)
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
double alphavalue(double gamma, TVector3 startpoint1, TVector3 startpoint2, TVector3 dircos1, TVector3 dircos2)
bool IsInVertexCollection(int a, std::vector< std::vector< int >> vertex_collection)
int IndexInVertexCollection(int a, int b, std::vector< std::vector< int >> vertex_collection)
1.8.11 
