From c7d002a056f11c0507a30b21326b31adcf1ba684 Mon Sep 17 00:00:00 2001 From: Valerio Pia Date: Tue, 17 Jun 2025 15:23:28 +0200 Subject: [PATCH 1/7] Added first vertexing version. It writes a vertices branch in the output --- CMakeLists.txt | 2 +- include/SANDTrackerVertexing.h | 7 +- include/StructLinkDef.h | 2 + include/struct.h | 9 ++ src/SANDDigitizationEDEPSIM.cpp | 8 +- src/SANDTrackerVertexing.cpp | 31 ++++- src/reconstruction.cpp | 194 +++++++++++++++++++++++++++----- tests/TestVertexing.cpp | 3 +- 8 files changed, 217 insertions(+), 39 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ebed57..dff0bc9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -91,7 +91,7 @@ target_link_libraries(Digitize Struct SANDGeoManager Utils) # Creates Reconstruct executable. add_executable(Reconstruct src/reconstruction.cpp) -target_link_libraries(Reconstruct PUBLIC Struct Utils SANDGeoManager SANDClustering EDEPTree TrackletFinder SANDTrackerCluster SANDTrackerDigit SANDTrackerUtils SANDKalmanFilter ROOT::EG) +target_link_libraries(Reconstruct PUBLIC Struct Utils SANDGeoManager SANDClustering EDEPTree TrackletFinder SANDTrackerCluster SANDTrackerDigit SANDTrackerUtils SANDKalmanFilter SANDTrackerVertexing ROOT::EG) # Creates Analyze executable. add_executable(Analyze src/analysis.cpp) diff --git a/include/SANDTrackerVertexing.h b/include/SANDTrackerVertexing.h index 9c70bdf..5ac30be 100644 --- a/include/SANDTrackerVertexing.h +++ b/include/SANDTrackerVertexing.h @@ -51,8 +51,11 @@ class Vertex { class TrackerVertexing { public: - void setParameters(double dz, double ip, double merging_radius, std::vector tracks); + void setParameters(double dz = 15, double ip = 15, double merging_radius = 15); + void setTracks(std::vector tracks); + void setTracks(std::vector tracks); int run(); + const std::vector& getVertices() {return vertices_;}; private: @@ -65,7 +68,7 @@ class TrackerVertexing { int mergeVertex(); void refineVertexPosition(double stepSize = 0.001, int nSteps = 50); - int id_vert_; + int id_vert_ = 0; double dz_; double ip_; double merging_radius_; diff --git a/include/StructLinkDef.h b/include/StructLinkDef.h index b45c078..ff42279 100644 --- a/include/StructLinkDef.h +++ b/include/StructLinkDef.h @@ -10,6 +10,7 @@ #pragma link C++ class std::vector < dg_ps> + ; #pragma link C++ class std::vector < dg_cell> + ; #pragma link C++ class std::vector < reco_cell > +; +#pragma link C++ class std::vector < vertex > +; //#pragma link C++ class std::map < std::string, std::vector < hit>> + ; #pragma link C++ class std::vector < dg_wire > +; #pragma link C++ class std::vector < track > +; @@ -27,4 +28,5 @@ #pragma link C++ class track + ; // #pragma link C++ class particle + ; #pragma link C++ class event + ; +#pragma link C++ class vertex + ; #endif diff --git a/include/struct.h b/include/struct.h index d119ef8..a9262e7 100644 --- a/include/struct.h +++ b/include/struct.h @@ -244,6 +244,15 @@ struct particle std::vector daughters; }; +struct vertex +{ + int id; + double x; + double y; + double z; + std::vector track_ids; +}; + struct event { double x; diff --git a/src/SANDDigitizationEDEPSIM.cpp b/src/SANDDigitizationEDEPSIM.cpp index 961ceb0..ee1f851 100644 --- a/src/SANDDigitizationEDEPSIM.cpp +++ b/src/SANDDigitizationEDEPSIM.cpp @@ -99,10 +99,10 @@ void CreateDigitsFromHits(const SANDGeoManager& geo, { long did = it->first(); // wire unique id const sand_geometry::tracker::WireInfo& wire_info = geo.getCellInfo(it->first())->second.getWire(); - double wire_time = 999.; - double drift_time = 999.; - double signal_time = 999.; - double t_hit = 999.; + double wire_time = 10e8; + double drift_time = 10e8; + double signal_time = 10e8; + double t_hit = 10e8; dg_wire d; d.det = it->second[0].det; diff --git a/src/SANDTrackerVertexing.cpp b/src/SANDTrackerVertexing.cpp index 3f4ecf4..7097b1d 100644 --- a/src/SANDTrackerVertexing.cpp +++ b/src/SANDTrackerVertexing.cpp @@ -321,11 +321,35 @@ void TrackerVertexing::flagVertex() { } } -void TrackerVertexing::setParameters(double dz, double ip, double merging_radius, std::vector tracks) { +// vertexing tracks +void TrackerVertexing::setTracks(std::vector tracks) { + clearAll(); + tracks_ = tracks; +} + +// struct tracks +void TrackerVertexing::setTracks(std::vector tracks) { + clearAll(); + for (const auto& t : tracks) { + Track converted_track; + converted_track.id = t.tid; + converted_track.x = t.x0; + converted_track.y = t.y0; + converted_track.z = t.z0; + converted_track.tx = t.b; + + double m_yz = (t.y0 - t.yc) / (t.z0 - t.zc); + converted_track.ty = m_yz; + std::cout << "Computed tan(phi) = " << converted_track.ty << std::endl; + + tracks_.push_back(converted_track); + } +} + +void TrackerVertexing::setParameters(double dz, double ip, double merging_radius) { dz_ = dz; ip_ = ip; merging_radius_ = merging_radius; - tracks_ = tracks; std::cout << tracks_.size() << " tracks read." << std::endl; std::cout << "DZ: " << dz_ << "\nIP: " << ip_ @@ -361,8 +385,7 @@ int TrackerVertexing::run() { std::cout << "Multi-Prong: " << vertices_multi_prong_.size() << std::endl; std::cout << "=========================" << std::endl; std::cout << "\nDump vertexes\n"; - dumpVertex("vertices.txt"); + // dumpVertex("vertices.txt"); - clearAll(); return 0; } diff --git a/src/reconstruction.cpp b/src/reconstruction.cpp index a57eab2..a133137 100644 --- a/src/reconstruction.cpp +++ b/src/reconstruction.cpp @@ -20,6 +20,7 @@ //added for kalman filter #include "SANDTrackletFinder.h" #include "SANDKalmanFilter.h" +#include "SANDTrackerVertexing.h" #include @@ -1604,12 +1605,91 @@ track runKalmanFilterManager(sand_reco::kf::TrackletMap z_to_tracklets, SParticl trk.b = reco_state.tanLambda(); trk.x0 = reco_state.x(); trk.y0 = reco_state.y(); - trk.z0 = reco_state.y(); + + trk.z0 = step.getZ(); + trk.yc = reco_state.y() - reco_state.radius() * sin(reco_state.phi()); + trk.zc = step.getZ() - reco_state.radius() * cos(reco_state.phi()); + std::cout << "tan(Phi) from KF = " << tan(reco_state.phi()) << std::endl; } return trk; } +void ProcessEventWithMC(std::vector& tracks, SANDGeoManager* sand_geo, TG4Event* mc_event, std::vector* digits, std::string tracker_name) +{ + + EDEPTree tree; + tree.InizializeFromEdep(*mc_event, sand_geo->getTGeoManager()); + + std::vector primaryTrj; + tree.Filter(std::back_insert_iterator>(primaryTrj), + [](const EDEPTrajectory& trj) { return trj.GetParentId() == -1;} ); + + TDatabasePDG pdg_db; + std::cout << "primaryTrj size: " << primaryTrj.size() << std::endl; + + for (auto trj:primaryTrj) { + + std::cout << "trj.GetPDGCode(): " << trj.GetPDGCode() << std::endl; + std::cout << "trj.GetId(): " << trj.GetId() << std::endl; + + auto particle = pdg_db.GetParticle(trj.GetPDGCode()); + + if (!particle) { + continue; + } + std::cout << "Found particle pdg" << std::endl; + + + if (particle->Mass() == 0 || particle->Charge() == 0) { + continue; + } + std::cout << "Particle in not neutral" << std::endl; + + // std::string log; + // trj.Print(log); + + for (auto h:trj.GetHitMap()) { + std::cout << component_to_string[h.first] << std::endl; + } + + if (trj.GetHitMap().find(string_to_component[tracker_name]) == trj.GetHitMap().end()) { + continue; + } + std::cout << "Found hits in tracker" << std::endl; + + + for (const auto& vertex:mc_event->Primaries) { + auto primary_trj_it = std::find_if(vertex.Particles.begin(), vertex.Particles.end(), [trj](TG4PrimaryParticle primary_trj){return primary_trj.GetTrackId() == trj.GetId();}); + if (primary_trj_it != vertex.Particles.end()) { + vertex.GetPosition().Print(); + break; + } + } + + + + auto trj_points = trj.GetTrajectoryPoints().at(string_to_component[tracker_name]); + auto state_vector = sand_reco::kf::utils::getStateVector(trj_points[1].GetMomentum(), + trj_points[1].GetPosition().Vect(), + particle->Charge()); + + track trk; + trk.tid = trj.GetId(); + trk.r = state_vector.radius(); + trk.h = state_vector.charge(); + trk.b = state_vector.tanLambda(); + trk.x0 = state_vector.x(); + trk.y0 = state_vector.y(); + + trk.z0 = trj_points[1].GetPosition().Z(); + trk.yc = state_vector.y() - state_vector.radius() * sin(state_vector.phi()); + trk.zc = trj_points[1].GetPosition().Z() - state_vector.radius() * cos(state_vector.phi()); + + tracks.push_back(trk); + } +} + void ProcessEventWithKF(std::vector& tracks, SANDGeoManager* sand_geo, TG4Event* mc_event, std::vector* digits) { int p[9] = {100, -2000, 2000, 100, -4000, -1000, 100, 23800, 26000}; @@ -1663,22 +1743,30 @@ void ProcessEventWithKF(std::vector& tracks, SANDGeoManager* sand_geo, TG [](const EDEPTrajectory& trj) { return trj.GetParentId() == -1;} ); TDatabasePDG pdg_db; + std::cout << "primaryTrj size: " << primaryTrj.size() << std::endl; + std::vector particleInfos; for (auto trj:primaryTrj) { - if (trj.GetHitMap().find(string_to_component[tracker_name]) == trj.GetHitMap().end()) { - continue; - } - + std::cout << "Particle pdg " << trj.GetPDGCode() << std::endl; auto particle = pdg_db.GetParticle(trj.GetPDGCode()); - if (!particle) { continue; } - + std::cout << "Found particle " << std::endl; + if (particle->Mass() == 0 || particle->Charge() == 0) { continue; } + std::cout << "Particle in not neutral/massless" << std::endl; + + for (auto h:trj.GetHitMap()) { + std::cout << component_to_string[h.first] << std::endl; + } + if (trj.GetHitMap().find(string_to_component[tracker_name]) == trj.GetHitMap().end()) { + continue; + } + std::cout << "Found hits in tracker" << std::endl; SParticleInfo pi; pi.pdg_code = trj.GetPDGCode(); @@ -1713,6 +1801,7 @@ void ProcessEventWithKF(std::vector& tracks, SANDGeoManager* sand_geo, TG } for (int ip = 0; ip < nParticles; ip++) { + std::cout << "nParticle: " << ip << std::endl; tracks.push_back(runKalmanFilterManager(z_to_tracklets, particleInfos[ip])); } } @@ -1721,7 +1810,8 @@ enum class STT_Mode { fast_only_primaries, fast, full, - primary_only_kf + primary_only_kf, + mc_primaries }; enum class ECAL_Mode { fast @@ -1729,7 +1819,7 @@ enum class ECAL_Mode { void Reconstruct(std::string const& fname_hits, std::string const& fname_digits, std::string const& fname_out, STT_Mode stt_mode, - ECAL_Mode ecal_mode) + ECAL_Mode ecal_mode, double dz, double ip, double mr) { std::cout << "Reconstruct\ninput hits: " << fname_hits << "\ninput digits: " << fname_digits @@ -1799,10 +1889,12 @@ void Reconstruct(std::string const& fname_hits, std::string const& fname_digits, std::vector vec_tr; std::vector vec_cl; + std::vector vec_vtx; TTree tout("tReco", "tReco"); tout.Branch("track", "std::vector", &vec_tr); tout.Branch("cluster", "std::vector", &vec_cl); + tout.Branch("vertex", "std::vector", &vec_vtx); const int nev = t->GetEntries(); const double epsilon = 0.5; @@ -1816,7 +1908,7 @@ void Reconstruct(std::string const& fname_hits, std::string const& fname_digits, std::cout << "Events: " << nev << " ["; std::cout << std::setw(3) << int(0) << "%]" << std::flush; - for (int i = 0; i < nev; i++) { + for (int i = 0; i < 1; i++) { t->GetEntry(i); std::cout << "\b\b\b\b\b" << std::setw(3) << int(double(i) / nev * 100) @@ -1851,8 +1943,30 @@ void Reconstruct(std::string const& fname_hits, std::string const& fname_digits, case STT_Mode::primary_only_kf: ProcessEventWithKF(vec_tr, &sand_geo, ev, vec_digi); break; + case STT_Mode::mc_primaries: + ProcessEventWithMC(vec_tr, &sand_geo, ev, vec_digi, trackerType); + break; + } + + TrackerVertexing vertex_finder; + vertex_finder.setTracks(vec_tr); + vertex_finder.setParameters(dz, ip, mr); + vertex_finder.run(); + auto vertices = vertex_finder.getVertices(); + + for (const auto& v:vertices) { + vertex vtx; + vtx.id = v.id; + vtx.x = v.x; + vtx.y = v.y; + vtx.z = v.z; + for (int j = 0; j < static_cast(v.indexTrack.size());j++) { + vtx.track_ids.push_back(vec_tr.at(v.indexTrack.at(j)).tid); + } + vec_vtx.push_back(vtx); } + switch (ecal_mode) { case ECAL_Mode::fast: // PreCluster(vec_cell, vec_cl); @@ -1867,6 +1981,7 @@ void Reconstruct(std::string const& fname_hits, std::string const& fname_digits, vec_tr.clear(); vec_cl.clear(); + vec_vtx.clear(); vec_digi->clear(); vec_cell->clear(); @@ -1882,35 +1997,60 @@ void Reconstruct(std::string const& fname_hits, std::string const& fname_digits, void help_reco() { std::cout - << "usage: Reconstruct hit_file digit_file output_file [stt_mode]\n"; + << "usage: Reconstruct hit_file digit_file output_file [stt_mode] [dz] [impact parameter] [merging_radius]\n"; std::cout << " - stt_mode: 'stt_mode::fast_only_primaries' (default) \n"; std::cout << " 'stt_mode::fast' \n"; std::cout << " 'stt_mode::full' \n"; std::cout << " 'stt_mode::primary_only_kf' \n"; + std::cout << " 'stt_mode::mc_primaries' \n"; + std::cout << " - dz: longitudinal track-vertex distance (default 15 mm)\n"; + std::cout << " - impact parameter: for 2-prongs vertices (default 15 mm)\n"; + std::cout << " - merging radius: radius for 2-prongs vertices clusterization (default 15 mm)\n"; } int main(int argc, char* argv[]) { // boost::program_options wuold be great here.... - if (argc < 4 || argc > 6) { + if (argc == 5) { + auto stt_mode = STT_Mode::fast_only_primaries; + if (argc > 4 && strcmp(argv[4], "stt_mode::full") == 0) { + stt_mode = STT_Mode::full; + std::cout << "STT_Mode: full\n"; + } else if (argc > 4 && strcmp(argv[4], "stt_mode::fast") == 0) { + stt_mode = STT_Mode::fast; + std::cout << "STT_Mode: fast\n"; + } else if (argc > 4 && strcmp(argv[4], "stt_mode::primary_only_kf") == 0) { + stt_mode = STT_Mode::primary_only_kf; + std::cout << "STT_Mode: kalman filter\n"; + } else if (argc > 4 && strcmp(argv[4], "stt_mode::mc_primaries") == 0) { + stt_mode = STT_Mode::mc_primaries; + std::cout << "STT_Mode: mc primaries\n"; + } else { + std::cout << "STT_Mode: fast_only_primaries\n"; + } + Reconstruct(argv[1], argv[2], argv[3], stt_mode, ECAL_Mode::fast, 15, 15, 15); + } else if (argc == 8 ) { + auto stt_mode = STT_Mode::fast_only_primaries; + if (argc > 4 && strcmp(argv[4], "stt_mode::full") == 0) { + stt_mode = STT_Mode::full; + std::cout << "STT_Mode: full\n"; + } else if (argc > 4 && strcmp(argv[4], "stt_mode::fast") == 0) { + stt_mode = STT_Mode::fast; + std::cout << "STT_Mode: fast\n"; + } else if (argc > 4 && strcmp(argv[4], "stt_mode::primary_only_kf") == 0) { + stt_mode = STT_Mode::primary_only_kf; + std::cout << "STT_Mode: kalman filter\n"; + } else if (argc > 4 && strcmp(argv[4], "stt_mode::mc_primaries") == 0) { + stt_mode = STT_Mode::mc_primaries; + std::cout << "STT_Mode: mc primaries\n"; + } else { + std::cout << "STT_Mode: fast_only_primaries\n"; + } + Reconstruct(argv[1], argv[2], argv[3], stt_mode, ECAL_Mode::fast, std::stod(argv[5]), std::stod(argv[6]), std::stod(argv[7])); + } else { help_reco(); return -1; } - - auto stt_mode = STT_Mode::fast_only_primaries; - if (argc > 4 && strcmp(argv[4], "stt_mode::full") == 0) { - stt_mode = STT_Mode::full; - std::cout << "STT_Mode: full\n"; - } else if (argc > 4 && strcmp(argv[4], "stt_mode::fast") == 0) { - stt_mode = STT_Mode::fast; - std::cout << "STT_Mode: fast\n"; - } else if (argc > 4 && strcmp(argv[4], "stt_mode::primary_only_kf") == 0) { - stt_mode = STT_Mode::primary_only_kf; - std::cout << "STT_Mode: kalman filter\n"; - } else { - std::cout << "STT_Mode: fast_only_primaries\n"; - } - Reconstruct(argv[1], argv[2], argv[3], stt_mode, ECAL_Mode::fast); return 0; } diff --git a/tests/TestVertexing.cpp b/tests/TestVertexing.cpp index 34696eb..cd79a04 100644 --- a/tests/TestVertexing.cpp +++ b/tests/TestVertexing.cpp @@ -70,7 +70,8 @@ int main(int argc, char* argv[]) std::cout << "Total tracks: " << tracks.size() << std::endl; TrackerVertexing vertex_finder; - vertex_finder.setParameters(std::stoi(argv[4]), std::stoi(argv[5]), std::stoi(argv[6]), tracks); + vertex_finder.setParameters(std::stoi(argv[4]), std::stoi(argv[5]), std::stoi(argv[6])); + vertex_finder.setTracks(tracks); vertex_finder.run(); } From 28cae9c5f059a4cea126e2efb7503f80f4376a93 Mon Sep 17 00:00:00 2001 From: Valerio Pia Date: Tue, 17 Jun 2025 15:25:11 +0200 Subject: [PATCH 2/7] Using first trajectory points --- src/reconstruction.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/reconstruction.cpp b/src/reconstruction.cpp index a133137..d715baf 100644 --- a/src/reconstruction.cpp +++ b/src/reconstruction.cpp @@ -1670,8 +1670,8 @@ void ProcessEventWithMC(std::vector& tracks, SANDGeoManager* sand_geo, TG auto trj_points = trj.GetTrajectoryPoints().at(string_to_component[tracker_name]); - auto state_vector = sand_reco::kf::utils::getStateVector(trj_points[1].GetMomentum(), - trj_points[1].GetPosition().Vect(), + auto state_vector = sand_reco::kf::utils::getStateVector(trj_points[0].GetMomentum(), + trj_points[0].GetPosition().Vect(), particle->Charge()); track trk; @@ -1682,9 +1682,9 @@ void ProcessEventWithMC(std::vector& tracks, SANDGeoManager* sand_geo, TG trk.x0 = state_vector.x(); trk.y0 = state_vector.y(); - trk.z0 = trj_points[1].GetPosition().Z(); + trk.z0 = trj_points[0].GetPosition().Z(); trk.yc = state_vector.y() - state_vector.radius() * sin(state_vector.phi()); - trk.zc = trj_points[1].GetPosition().Z() - state_vector.radius() * cos(state_vector.phi()); + trk.zc = trj_points[0].GetPosition().Z() - state_vector.radius() * cos(state_vector.phi()); tracks.push_back(trk); } From b20e1639bc9084d152323e253167fd7255bca88a Mon Sep 17 00:00:00 2001 From: Valerio Pia Date: Tue, 17 Jun 2025 15:35:10 +0200 Subject: [PATCH 3/7] Fixed misseed clear, processing all events --- src/reconstruction.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/reconstruction.cpp b/src/reconstruction.cpp index d715baf..544cced 100644 --- a/src/reconstruction.cpp +++ b/src/reconstruction.cpp @@ -1908,7 +1908,7 @@ void Reconstruct(std::string const& fname_hits, std::string const& fname_digits, std::cout << "Events: " << nev << " ["; std::cout << std::setw(3) << int(0) << "%]" << std::flush; - for (int i = 0; i < 1; i++) { + for (int i = 0; i < nev; i++) { t->GetEntry(i); std::cout << "\b\b\b\b\b" << std::setw(3) << int(double(i) / nev * 100) @@ -1916,6 +1916,7 @@ void Reconstruct(std::string const& fname_hits, std::string const& fname_digits, vec_tr.clear(); vec_cl.clear(); + vec_vtx.clear(); double xvtx_reco, yvtx_reco, zvtx_reco; int VtxType; From 4e8bb933f1de1ac9779e1ae903ad120e40f0ea88 Mon Sep 17 00:00:00 2001 From: Valerio Pia Date: Tue, 17 Jun 2025 16:26:25 +0200 Subject: [PATCH 4/7] Fixed bug with hor variable in dg_wire --- include/SANDTrackerPlane.h | 7 ++++++- src/SANDDigitizationEDEPSIM.cpp | 5 ++++- src/SANDGeoManager.cpp | 9 +++++++-- 3 files changed, 17 insertions(+), 4 deletions(-) diff --git a/include/SANDTrackerPlane.h b/include/SANDTrackerPlane.h index 7aad9fa..105bfac 100644 --- a/include/SANDTrackerPlane.h +++ b/include/SANDTrackerPlane.h @@ -74,9 +74,14 @@ class Plane void addCell(const double transverse_coordinate, Cell c) { if(coord_to_id_map_.find(transverse_coordinate) == coord_to_id_map_.end()) { - c.setPlane(this); coord_to_id_map_.insert({transverse_coordinate, c.getId()}); id_to_cell_map_.insert({c.getId(), c}); + id_to_cell_map_[c.getId()].setPlane(this); + } + } + void updateCells() { + for (auto& c:id_to_cell_map_) { + c.second.setPlane(this); } } std::map& getIdToCellMap() {return id_to_cell_map_;}; diff --git a/src/SANDDigitizationEDEPSIM.cpp b/src/SANDDigitizationEDEPSIM.cpp index ee1f851..02ce110 100644 --- a/src/SANDDigitizationEDEPSIM.cpp +++ b/src/SANDDigitizationEDEPSIM.cpp @@ -98,7 +98,8 @@ void CreateDigitsFromHits(const SANDGeoManager& geo, it != hits2cell.end(); ++it) // run over wires { long did = it->first(); // wire unique id - const sand_geometry::tracker::WireInfo& wire_info = geo.getCellInfo(it->first())->second.getWire(); + const auto& cell_info = geo.getCellInfo(it->first())->second; + const sand_geometry::tracker::WireInfo& wire_info = cell_info.getWire(); double wire_time = 10e8; double drift_time = 10e8; double signal_time = 10e8; @@ -107,6 +108,8 @@ void CreateDigitsFromHits(const SANDGeoManager& geo, dg_wire d; d.det = it->second[0].det; d.did = did; + if(!cell_info.getPlane()) std::cout << "Sto male" << std::endl; + d.hor = (cell_info.getPlane()->getRotation() == 0) ? true : false; d.de = 0; // To Do: what point do we want to save? // Center or one of the attachment points? diff --git a/src/SANDGeoManager.cpp b/src/SANDGeoManager.cpp index 2668a24..e659c98 100644 --- a/src/SANDGeoManager.cpp +++ b/src/SANDGeoManager.cpp @@ -1519,8 +1519,13 @@ void SANDGeoManager::rearrangePlanes() for (auto it = planes_.begin(); it != planes_.end(); ++it) { - id_to_plane_[it->uId()] = it; - } + id_to_plane_[it->uId()] = it; + } + + for (auto it = planes_.begin(); + it != planes_.end(); ++it) { + it->updateCells(); + } } void SANDGeoManager::setTrackerInfo() From 07870e895e2195c8072c3f8c2023ac85f37b5cb4 Mon Sep 17 00:00:00 2001 From: Valerio Pia Date: Tue, 8 Jul 2025 11:02:33 +0200 Subject: [PATCH 5/7] Fixed CMake --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index dff0bc9..efeded6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -153,7 +153,7 @@ endif() # Copy setup.sh configuration file configure_file(setup.sh "${CMAKE_INSTALL_PREFIX}/setup.sh" COPYONLY) - install(TARGETS Utils Struct SANDEventDisplay SANDGeoManager Digitize Reconstruct Display FastCheck eventDisplay TestKF TestSeed TestVertexing SANDClustering SANDECALClustering NoiseMatrixComputation TrackletFinder SANDTrackerCluster SANDTrackerDigit SANDKalmanFilter SANDKFTrack SANDTrackerUtils + install(TARGETS Utils Struct SANDEventDisplay SANDGeoManager Digitize Reconstruct Display FastCheck eventDisplay TestKF TestSeed TestVertexing SANDClustering SANDECALClustering NoiseMatrixComputation TrackletFinder SANDTrackerCluster SANDTrackerDigit SANDKalmanFilter SANDKFTrack SANDTrackerUtils SANDTrackerVertexing EXPORT SandRecoTargets RUNTIME DESTINATION "${CMAKE_INSTALL_PREFIX}/bin" LIBRARY DESTINATION "${CMAKE_INSTALL_PREFIX}/lib" From 1f1825d7073e59d5bb244834af104f8b5c13bd8b Mon Sep 17 00:00:00 2001 From: Valerio Pia Date: Thu, 17 Jul 2025 12:27:11 +0200 Subject: [PATCH 6/7] Removed couts, fixed argument parsing --- src/SANDTrackerVertexing.cpp | 29 +++++++------ src/reconstruction.cpp | 80 ++++++++++++++++++------------------ 2 files changed, 54 insertions(+), 55 deletions(-) diff --git a/src/SANDTrackerVertexing.cpp b/src/SANDTrackerVertexing.cpp index 7097b1d..6c7d981 100644 --- a/src/SANDTrackerVertexing.cpp +++ b/src/SANDTrackerVertexing.cpp @@ -180,7 +180,7 @@ double pointTo3DlineDistance(double x0, double y0, double z0, Track tr) { int TrackerVertexing::mergeVertex() { int nVertex = static_cast(vertices_list_.size()); if (nVertex == 0) { - std::cout << "No vertex to merge" << std::endl; + // std::cout << "No vertex to merge" << std::endl; return 0; } for (int i = 0; i < nVertex; i++) { @@ -340,7 +340,6 @@ void TrackerVertexing::setTracks(std::vector tracks) { double m_yz = (t.y0 - t.yc) / (t.z0 - t.zc); converted_track.ty = m_yz; - std::cout << "Computed tan(phi) = " << converted_track.ty << std::endl; tracks_.push_back(converted_track); } @@ -351,19 +350,19 @@ void TrackerVertexing::setParameters(double dz, double ip, double merging_radius ip_ = ip; merging_radius_ = merging_radius; - std::cout << tracks_.size() << " tracks read." << std::endl; - std::cout << "DZ: " << dz_ << "\nIP: " << ip_ - << "\nMerging Radius: " << merging_radius << std::endl; + // std::cout << tracks_.size() << " tracks read." << std::endl; + // std::cout << "DZ: " << dz_ << "\nIP: " << ip_ + // << "\nMerging Radius: " << merging_radius << std::endl; } int TrackerVertexing::run() { int pairVertexes = doVertex(); - dumpVertex("2-Prong.txt"); - std::cout << "Start selecting neighboor vertexes (" << pairVertexes - << " - 2 Prong )" << std::endl; + // dumpVertex("2-Prong.txt"); + // std::cout << "Start selecting neighboor vertexes (" << pairVertexes + // << " - 2 Prong )" << std::endl; while (static_cast(vertices_.size()) > 0) selectVertex(); - std::cout << "Start merging..." << std::endl; + // std::cout << "Start merging..." << std::endl; mergeVertex(); refineVertexPosition(); @@ -377,14 +376,14 @@ int TrackerVertexing::run() { vertices_.push_back(vertices_2_prong_.at(j)); } - std::cout << "Flag vertexes\n"; + // std::cout << "Flag vertexes\n"; flagVertex(); - std::cout << "\n\n======== Results ========" << std::endl; - std::cout << "2-Prong: " << vertices_2_prong_.size() << std::endl; - std::cout << "Multi-Prong: " << vertices_multi_prong_.size() << std::endl; - std::cout << "=========================" << std::endl; - std::cout << "\nDump vertexes\n"; + // std::cout << "\n\n======== Results ========" << std::endl; + // std::cout << "2-Prong: " << vertices_2_prong_.size() << std::endl; + // std::cout << "Multi-Prong: " << vertices_multi_prong_.size() << std::endl; + // std::cout << "=========================" << std::endl; + // std::cout << "\nDump vertexes\n"; // dumpVertex("vertices.txt"); return 0; diff --git a/src/reconstruction.cpp b/src/reconstruction.cpp index 544cced..1ec6acb 100644 --- a/src/reconstruction.cpp +++ b/src/reconstruction.cpp @@ -1998,60 +1998,60 @@ void Reconstruct(std::string const& fname_hits, std::string const& fname_digits, void help_reco() { std::cout - << "usage: Reconstruct hit_file digit_file output_file [stt_mode] [dz] [impact parameter] [merging_radius]\n"; + << "usage: Reconstruct hit_file digit_file output_file [stt_mode] dz [dz value] impact_parameter [impact_parameter value] merging_radius [merging_radius value]\n"; std::cout << " - stt_mode: 'stt_mode::fast_only_primaries' (default) \n"; std::cout << " 'stt_mode::fast' \n"; std::cout << " 'stt_mode::full' \n"; std::cout << " 'stt_mode::primary_only_kf' \n"; std::cout << " 'stt_mode::mc_primaries' \n"; std::cout << " - dz: longitudinal track-vertex distance (default 15 mm)\n"; - std::cout << " - impact parameter: for 2-prongs vertices (default 15 mm)\n"; - std::cout << " - merging radius: radius for 2-prongs vertices clusterization (default 15 mm)\n"; + std::cout << " - impact_parameter: for 2-prongs vertices (default 15 mm)\n"; + std::cout << " - merging_radius: radius for 2-prongs vertices clusterization (default 15 mm)\n"; } int main(int argc, char* argv[]) { + if (argc < 4) { + help_reco(); + return -1; + } + // boost::program_options wuold be great here.... + auto stt_mode = STT_Mode::fast_only_primaries; + double dz = 15; + double impact_parameter = 15; + double merging_radius = 15; + + for (int i = 0; i < argc; i++) { + if (strcmp(argv[i], "stt_mode") == 0) { + if (strcmp(argv[i], "stt_mode::full") == 0) { + stt_mode = STT_Mode::full; + std::cout << "STT_Mode: full\n"; + } else if (strcmp(argv[i], "stt_mode::fast") == 0) { + stt_mode = STT_Mode::fast; + std::cout << "STT_Mode: fast\n"; + } else if (strcmp(argv[i], "stt_mode::primary_only_kf") == 0) { + stt_mode = STT_Mode::primary_only_kf; + std::cout << "STT_Mode: kalman filter\n"; + } else if (strcmp(argv[i], "stt_mode::mc_primaries") == 0) { + stt_mode = STT_Mode::mc_primaries; + std::cout << "STT_Mode: mc primaries\n"; + } else { + std::cout << "STT_Mode: fast_only_primaries\n"; + } + } - if (argc == 5) { - auto stt_mode = STT_Mode::fast_only_primaries; - if (argc > 4 && strcmp(argv[4], "stt_mode::full") == 0) { - stt_mode = STT_Mode::full; - std::cout << "STT_Mode: full\n"; - } else if (argc > 4 && strcmp(argv[4], "stt_mode::fast") == 0) { - stt_mode = STT_Mode::fast; - std::cout << "STT_Mode: fast\n"; - } else if (argc > 4 && strcmp(argv[4], "stt_mode::primary_only_kf") == 0) { - stt_mode = STT_Mode::primary_only_kf; - std::cout << "STT_Mode: kalman filter\n"; - } else if (argc > 4 && strcmp(argv[4], "stt_mode::mc_primaries") == 0) { - stt_mode = STT_Mode::mc_primaries; - std::cout << "STT_Mode: mc primaries\n"; - } else { - std::cout << "STT_Mode: fast_only_primaries\n"; + if (strcmp(argv[i], "dz") == 0) { + dz = std::stod(argv[i+1]); } - Reconstruct(argv[1], argv[2], argv[3], stt_mode, ECAL_Mode::fast, 15, 15, 15); - } else if (argc == 8 ) { - auto stt_mode = STT_Mode::fast_only_primaries; - if (argc > 4 && strcmp(argv[4], "stt_mode::full") == 0) { - stt_mode = STT_Mode::full; - std::cout << "STT_Mode: full\n"; - } else if (argc > 4 && strcmp(argv[4], "stt_mode::fast") == 0) { - stt_mode = STT_Mode::fast; - std::cout << "STT_Mode: fast\n"; - } else if (argc > 4 && strcmp(argv[4], "stt_mode::primary_only_kf") == 0) { - stt_mode = STT_Mode::primary_only_kf; - std::cout << "STT_Mode: kalman filter\n"; - } else if (argc > 4 && strcmp(argv[4], "stt_mode::mc_primaries") == 0) { - stt_mode = STT_Mode::mc_primaries; - std::cout << "STT_Mode: mc primaries\n"; - } else { - std::cout << "STT_Mode: fast_only_primaries\n"; + if (strcmp(argv[i], "impact_parameter") == 0) { + impact_parameter = std::stod(argv[i+1]); + } + if (strcmp(argv[i], "merging_radius") == 0) { + merging_radius = std::stod(argv[i+1]); } - Reconstruct(argv[1], argv[2], argv[3], stt_mode, ECAL_Mode::fast, std::stod(argv[5]), std::stod(argv[6]), std::stod(argv[7])); - } else { - help_reco(); - return -1; } + + Reconstruct(argv[1], argv[2], argv[3], stt_mode, ECAL_Mode::fast, dz, impact_parameter, merging_radius); return 0; } From 7c78948f68e15f2db70d53947a0c4b159b5a11c1 Mon Sep 17 00:00:00 2001 From: Valerio Pia Date: Tue, 16 Sep 2025 15:13:25 +0200 Subject: [PATCH 7/7] Removed unneded couts, changed some constant values to std::max --- src/SANDDigitizationEDEPSIM.cpp | 14 +++++++++----- src/reconstruction.cpp | 28 ---------------------------- 2 files changed, 9 insertions(+), 33 deletions(-) diff --git a/src/SANDDigitizationEDEPSIM.cpp b/src/SANDDigitizationEDEPSIM.cpp index 02ce110..2df4ce2 100644 --- a/src/SANDDigitizationEDEPSIM.cpp +++ b/src/SANDDigitizationEDEPSIM.cpp @@ -3,6 +3,7 @@ #include #include +#include #include @@ -100,15 +101,18 @@ void CreateDigitsFromHits(const SANDGeoManager& geo, long did = it->first(); // wire unique id const auto& cell_info = geo.getCellInfo(it->first())->second; const sand_geometry::tracker::WireInfo& wire_info = cell_info.getWire(); - double wire_time = 10e8; - double drift_time = 10e8; - double signal_time = 10e8; - double t_hit = 10e8; + double wire_time = std::numeric_limits::max(); + double drift_time = std::numeric_limits::max(); + double signal_time = std::numeric_limits::max(); + double t_hit = std::numeric_limits::max(); dg_wire d; d.det = it->second[0].det; d.did = did; - if(!cell_info.getPlane()) std::cout << "Sto male" << std::endl; + if(!cell_info.getPlane()) { + std::cerr << "Plane not found. Skipping hit.." << std::endl; + continue; + } d.hor = (cell_info.getPlane()->getRotation() == 0) ? true : false; d.de = 0; // To Do: what point do we want to save? diff --git a/src/reconstruction.cpp b/src/reconstruction.cpp index 1ec6acb..5f2f4a3 100644 --- a/src/reconstruction.cpp +++ b/src/reconstruction.cpp @@ -1626,38 +1626,22 @@ void ProcessEventWithMC(std::vector& tracks, SANDGeoManager* sand_geo, TG [](const EDEPTrajectory& trj) { return trj.GetParentId() == -1;} ); TDatabasePDG pdg_db; - std::cout << "primaryTrj size: " << primaryTrj.size() << std::endl; for (auto trj:primaryTrj) { - std::cout << "trj.GetPDGCode(): " << trj.GetPDGCode() << std::endl; - std::cout << "trj.GetId(): " << trj.GetId() << std::endl; - auto particle = pdg_db.GetParticle(trj.GetPDGCode()); if (!particle) { continue; } - std::cout << "Found particle pdg" << std::endl; - if (particle->Mass() == 0 || particle->Charge() == 0) { continue; } - std::cout << "Particle in not neutral" << std::endl; - - // std::string log; - // trj.Print(log); - - for (auto h:trj.GetHitMap()) { - std::cout << component_to_string[h.first] << std::endl; - } if (trj.GetHitMap().find(string_to_component[tracker_name]) == trj.GetHitMap().end()) { continue; } - std::cout << "Found hits in tracker" << std::endl; - for (const auto& vertex:mc_event->Primaries) { auto primary_trj_it = std::find_if(vertex.Particles.begin(), vertex.Particles.end(), [trj](TG4PrimaryParticle primary_trj){return primary_trj.GetTrackId() == trj.GetId();}); @@ -1666,8 +1650,6 @@ void ProcessEventWithMC(std::vector& tracks, SANDGeoManager* sand_geo, TG break; } } - - auto trj_points = trj.GetTrajectoryPoints().at(string_to_component[tracker_name]); auto state_vector = sand_reco::kf::utils::getStateVector(trj_points[0].GetMomentum(), @@ -1743,30 +1725,21 @@ void ProcessEventWithKF(std::vector& tracks, SANDGeoManager* sand_geo, TG [](const EDEPTrajectory& trj) { return trj.GetParentId() == -1;} ); TDatabasePDG pdg_db; - std::cout << "primaryTrj size: " << primaryTrj.size() << std::endl; std::vector particleInfos; for (auto trj:primaryTrj) { - - std::cout << "Particle pdg " << trj.GetPDGCode() << std::endl; auto particle = pdg_db.GetParticle(trj.GetPDGCode()); if (!particle) { continue; } - std::cout << "Found particle " << std::endl; if (particle->Mass() == 0 || particle->Charge() == 0) { continue; } - std::cout << "Particle in not neutral/massless" << std::endl; - for (auto h:trj.GetHitMap()) { - std::cout << component_to_string[h.first] << std::endl; - } if (trj.GetHitMap().find(string_to_component[tracker_name]) == trj.GetHitMap().end()) { continue; } - std::cout << "Found hits in tracker" << std::endl; SParticleInfo pi; pi.pdg_code = trj.GetPDGCode(); @@ -1801,7 +1774,6 @@ void ProcessEventWithKF(std::vector& tracks, SANDGeoManager* sand_geo, TG } for (int ip = 0; ip < nParticles; ip++) { - std::cout << "nParticle: " << ip << std::endl; tracks.push_back(runKalmanFilterManager(z_to_tracklets, particleInfos[ip])); } }