diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ebed57..efeded6 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) @@ -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" 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/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..2df4ce2 100644 --- a/src/SANDDigitizationEDEPSIM.cpp +++ b/src/SANDDigitizationEDEPSIM.cpp @@ -3,6 +3,7 @@ #include #include +#include #include @@ -98,15 +99,21 @@ 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(); - double wire_time = 999.; - double drift_time = 999.; - double signal_time = 999.; - double t_hit = 999.; + const auto& cell_info = geo.getCellInfo(it->first())->second; + const sand_geometry::tracker::WireInfo& wire_info = cell_info.getWire(); + 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::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? // Center or one of the attachment points? diff --git a/src/SANDGeoManager.cpp b/src/SANDGeoManager.cpp index 29c0802..6fcb552 100644 --- a/src/SANDGeoManager.cpp +++ b/src/SANDGeoManager.cpp @@ -1516,8 +1516,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() diff --git a/src/SANDTrackerVertexing.cpp b/src/SANDTrackerVertexing.cpp index 3f4ecf4..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++) { @@ -321,25 +321,48 @@ 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; + + 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_ - << "\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(); @@ -353,16 +376,15 @@ 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"; - dumpVertex("vertices.txt"); + // 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"); - clearAll(); return 0; } diff --git a/src/reconstruction.cpp b/src/reconstruction.cpp index a57eab2..5f2f4a3 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,73 @@ 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; + + for (auto trj:primaryTrj) { + + auto particle = pdg_db.GetParticle(trj.GetPDGCode()); + + if (!particle) { + continue; + } + + if (particle->Mass() == 0 || particle->Charge() == 0) { + continue; + } + + if (trj.GetHitMap().find(string_to_component[tracker_name]) == trj.GetHitMap().end()) { + continue; + } + + 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[0].GetMomentum(), + trj_points[0].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[0].GetPosition().Z(); + trk.yc = state_vector.y() - state_vector.radius() * sin(state_vector.phi()); + trk.zc = trj_points[0].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,23 +1725,22 @@ void ProcessEventWithKF(std::vector& tracks, SANDGeoManager* sand_geo, TG [](const EDEPTrajectory& trj) { return trj.GetParentId() == -1;} ); TDatabasePDG pdg_db; + std::vector particleInfos; for (auto trj:primaryTrj) { - - if (trj.GetHitMap().find(string_to_component[tracker_name]) == trj.GetHitMap().end()) { - continue; - } - auto particle = pdg_db.GetParticle(trj.GetPDGCode()); - if (!particle) { continue; } - + if (particle->Mass() == 0 || particle->Charge() == 0) { continue; } + if (trj.GetHitMap().find(string_to_component[tracker_name]) == trj.GetHitMap().end()) { + continue; + } + SParticleInfo pi; pi.pdg_code = trj.GetPDGCode(); pi.id = trj.GetId(); @@ -1721,7 +1782,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 +1791,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 +1861,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; @@ -1824,6 +1888,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; @@ -1851,7 +1916,29 @@ 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: @@ -1867,6 +1954,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 +1970,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 [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"; } int main(int argc, char* argv[]) { - // boost::program_options wuold be great here.... - - if (argc < 4 || argc > 6) { + if (argc < 4) { help_reco(); return -1; } + // boost::program_options wuold be great here.... 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"; + 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 (strcmp(argv[i], "dz") == 0) { + dz = std::stod(argv[i+1]); + } + 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); + + Reconstruct(argv[1], argv[2], argv[3], stt_mode, ECAL_Mode::fast, dz, impact_parameter, merging_radius); 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(); }