From 48eb8dad8707ae22ab47c1ba289e81ceede1ea26 Mon Sep 17 00:00:00 2001
From: Riccardo DAmico <rdamico@neutrino-01.novalocal>
Date: Thu, 15 Feb 2024 11:06:57 +0100
Subject: [PATCH 01/12] Added GetFirstPointInDetector and
GetLastPointInDetector methods and first_points and last_points private
variables to EDEPTrajectory
---
include/EDEPTrajectory.h | 6 ++
src/EDEPTrajectory.cpp | 170 ++++++++++++++++++++++++++++++---------
2 files changed, 139 insertions(+), 37 deletions(-)
diff --git a/include/EDEPTrajectory.h b/include/EDEPTrajectory.h
index 62f8291..7d8d10c 100644
--- a/include/EDEPTrajectory.h
+++ b/include/EDEPTrajectory.h
@@ -42,6 +42,10 @@ class EDEPTrajectory {
int GetDepth() const {return depth_;};
int GetParentId() const {return parent_id_;};
int GetPDGCode() const {return pdg_code_;};
+
+ EDEPTrajectoryHits& GetFirstPointInDetector(component component_name);
+ EDEPTrajectoryHits& GetLastPointInDetector(component component_name);
+
TLorentzVector GetInitialMomentum() const {return p0_;};
std::vector<EDEPTrajectory>& GetChildrenTrajectories() {return children_trajectories_;};
@@ -102,6 +106,8 @@ class EDEPTrajectory {
std::vector<EDEPTrajectory> children_trajectories_;
std::map<component, bool> exiting_map_;
std::map<component, bool> entering_map_;
+ std::map<component, EDEPTrajectoryHits> last_points_;
+ std::map<component, EDEPTrajectoryHits> first_points_;
EDEPTrajectory* parent_trajectory_;
int id_;
int parent_id_;
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index fdc4782..5034439 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -28,42 +28,115 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
pdg_code_(trajectory.GetPDGCode()) {
for (auto it = trajectory.Points.begin(); it != trajectory.Points.end(); ++it) {
- trajectory_points_.push_back(*it);
-
- auto next_it = std::next(it);
-
- if (next_it == trajectory.Points.end()) {
- continue;
- }
-
- auto current_node = GetNode(*it, geo);
- auto current_volume = current_node->GetName();
- auto next_node = GetNode(*next_it, geo);
- auto next_volume = next_node->GetName();
-
- bool in_grain = Match(current_volume, grain_names);
- bool in_stt = Match(current_volume, stt_names);
- bool in_ecal = Match(current_volume, ecal_names);
- bool in_mag = Match(current_volume, magnet_names);
- bool in_world = Match(current_volume, world_names);
-
- bool next_grain = Match(next_volume, grain_names);
- bool next_stt = Match(next_volume, stt_names);
- bool next_ecal = Match(next_volume, ecal_names);
- bool next_mag = Match(next_volume, magnet_names);
- bool next_world = Match(next_volume, world_names);
-
- if (!exiting_map_[component::GRAIN]) exiting_map_[component::GRAIN] = (in_grain && ( next_stt || next_ecal || next_mag || next_world ));
- if (!exiting_map_[component::STRAW]) exiting_map_[component::STRAW] = (in_stt && ( next_grain || next_ecal || next_mag || next_world ));
- if (!exiting_map_[component::ECAL]) exiting_map_[component::ECAL] = (in_ecal && ( next_stt || next_grain || next_mag || next_world ));
- if (!exiting_map_[component::MAGNET]) exiting_map_[component::MAGNET] = (in_mag && ( next_stt || next_ecal || next_grain || next_world ));
- if (!exiting_map_[component::WORLD]) exiting_map_[component::WORLD] = (in_world && ( next_stt || next_ecal || next_mag || next_grain ));
-
- if (!entering_map_[component::GRAIN]) entering_map_[component::GRAIN] = (next_grain && ( in_stt || in_ecal || in_mag || in_world ));
- if (!entering_map_[component::STRAW]) entering_map_[component::STRAW] = (next_stt && ( in_grain || in_ecal || in_mag || in_world ));
- if (!entering_map_[component::ECAL]) entering_map_[component::ECAL] = (next_ecal && ( in_stt || in_grain || in_mag || in_world ));
- if (!entering_map_[component::MAGNET]) entering_map_[component::MAGNET] = (next_mag && ( in_stt || in_ecal || in_grain || in_world ));
- if (!entering_map_[component::WORLD]) entering_map_[component::WORLD] = (next_world && ( in_stt || in_ecal || in_mag || in_grain ));
+ trajectory_points_.push_back(*it);
+
+ auto next_it = std::next(it);
+
+ if (next_it == trajectory.Points.end()) {
+ continue;
+ }
+
+ auto current_node = GetNode(*it, geo);
+ auto current_volume = current_node->GetName();
+ auto next_node = GetNode(*next_it, geo);
+ auto next_volume = next_node->GetName();
+
+ bool in_grain = Match(current_volume, grain_names);
+ bool in_stt = Match(current_volume, stt_names);
+ bool in_ecal = Match(current_volume, ecal_names);
+ bool in_mag = Match(current_volume, magnet_names);
+ bool in_world = Match(current_volume, world_names);
+
+ bool next_grain = Match(next_volume, grain_names);
+ bool next_stt = Match(next_volume, stt_names);
+ bool next_ecal = Match(next_volume, ecal_names);
+ bool next_mag = Match(next_volume, magnet_names);
+ bool next_world = Match(next_volume, world_names);
+
+ if (!exiting_map_[component::GRAIN]) exiting_map_[component::GRAIN] = (in_grain && ( next_stt || next_ecal || next_mag || next_world ));
+ if (!exiting_map_[component::STRAW]) exiting_map_[component::STRAW] = (in_stt && ( next_grain || next_ecal || next_mag || next_world ));
+ if (!exiting_map_[component::ECAL]) exiting_map_[component::ECAL] = (in_ecal && ( next_stt || next_grain || next_mag || next_world ));
+ if (!exiting_map_[component::MAGNET]) exiting_map_[component::MAGNET] = (in_mag && ( next_stt || next_ecal || next_grain || next_world ));
+ if (!exiting_map_[component::WORLD]) exiting_map_[component::WORLD] = (in_world && ( next_stt || next_ecal || next_mag || next_grain ));
+
+ if (!entering_map_[component::GRAIN]) entering_map_[component::GRAIN] = (next_grain && ( in_stt || in_ecal || in_mag || in_world ));
+ if (!entering_map_[component::STRAW]) entering_map_[component::STRAW] = (next_stt && ( in_grain || in_ecal || in_mag || in_world ));
+ if (!entering_map_[component::ECAL]) entering_map_[component::ECAL] = (next_ecal && ( in_stt || in_grain || in_mag || in_world ));
+ if (!entering_map_[component::MAGNET]) entering_map_[component::MAGNET] = (next_mag && ( in_stt || in_ecal || in_grain || in_world ));
+ if (!entering_map_[component::WORLD]) entering_map_[component::WORLD] = (next_world && ( in_stt || in_ecal || in_mag || in_grain ));
+
+
+ if(in_grain && next_stt){
+ last_points_[component::GRAIN].push_back(*it);
+ first_points_[component::STRAW].push_back(*next_it);
+ } else if(in_grain && next_ecal){
+ last_points_[component::GRAIN].push_back(*it);
+ first_points_[component::ECAL].push_back(*next_it);
+ } else if(in_grain && next_mag){
+ last_points_[component::GRAIN].push_back(*it);
+ first_points_[component::MAGNET].push_back(*next_it);
+ } else if(in_grain && next_world){
+ last_points_[component::GRAIN].push_back(*it);
+ first_points_[component::WORLD].push_back(*next_it);
+ }
+
+ if(in_stt && next_grain){
+ last_points_[component::STRAW].push_back(*it);
+ first_points_[component::GRAIN].push_back(*next_it);
+ } else if(in_stt && next_ecal){
+ last_points_[component::STRAW].push_back(*it);
+ first_points_[component::ECAL].push_back(*next_it);
+ } else if(in_stt && next_mag){
+ last_points_[component::STRAW].push_back(*it);
+ first_points_[component::MAGNET].push_back(*next_it);
+ } else if(in_stt && next_world){
+ last_points_[component::STRAW].push_back(*it);
+ first_points_[component::WORLD].push_back(*next_it);
+ }
+
+ if(in_ecal && next_grain){
+ last_points_[component::ECAL].push_back(*it);
+ first_points_[component::GRAIN].push_back(*next_it);
+ } else if(in_ecal && next_stt){
+ last_points_[component::ECAL].push_back(*it);
+ first_points_[component::STRAW].push_back(*next_it);
+
+ } else if(in_ecal && next_mag){
+ last_points_[component::ECAL].push_back(*it);
+ first_points_[component::MAGNET].push_back(*next_it);
+
+ } else if(in_ecal && next_world){
+ last_points_[component::ECAL].push_back(*it);
+ first_points_[component::WORLD].push_back(*next_it);
+ }
+
+ if(in_mag && next_grain){
+ last_points_[component::MAGNET].push_back(*it);
+ first_points_[component::GRAIN].push_back(*next_it);
+ } else if(in_mag && next_stt){
+ last_points_[component::MAGNET].push_back(*it);
+ first_points_[component::STRAW].push_back(*next_it);
+ } else if(in_mag && next_ecal){
+ last_points_[component::MAGNET].push_back(*it);
+ first_points_[component::ECAL].push_back(*next_it);
+ } else if(in_mag && next_world){
+ last_points_[component::MAGNET].push_back(*it);
+ first_points_[component::WORLD].push_back(*next_it);
+ }
+
+ if(in_world && next_grain){
+ last_points_[component::WORLD].push_back(*it);
+ first_points_[component::GRAIN].push_back(*next_it);
+ } else if(in_world && next_stt){
+ last_points_[component::WORLD].push_back(*it);
+ first_points_[component::STRAW].push_back(*next_it);
+ } else if(in_world && next_ecal){
+ last_points_[component::WORLD].push_back(*it);
+ first_points_[component::ECAL].push_back(*next_it);
+ } else if(in_world && next_world){
+ last_points_[component::WORLD].push_back(*it);
+ first_points_[component::MAGNET].push_back(*next_it);
+ }
}
for (const auto& hmap:hits) {
@@ -80,6 +153,8 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
hit_map_[string_to_component[hmap.first]] = edep_hits;
}
}
+
+
}
EDEPTrajectory::EDEPTrajectory(const EDEPTrajectory& trj) {
@@ -100,7 +175,9 @@ EDEPTrajectory& EDEPTrajectory::operator=(const EDEPTrajectory& trj) {
this->hit_map_ = trj.hit_map_;
this->entering_map_ = trj.entering_map_;
this->exiting_map_ = trj.exiting_map_;
- this->trajectory_points_ = trj.trajectory_points_;
+ this -> last_points_ = trj.last_points_;
+ this-> first_points_ = trj.first_points_;
+ this->trajectory_points_ = trj.trajectory_points_;
return *this;
}
@@ -115,6 +192,8 @@ bool EDEPTrajectory::operator==(const EDEPTrajectory& trj) {
// this->hit_map_ == trj.hit_map_ &&
this->entering_map_ == trj.entering_map_ &&
this->exiting_map_ == trj.exiting_map_
+ // this->last_points_ == trj.last_points_ &&
+ // this -> first_points_ == trj.first_points_
// this->trajectory_points_ == trj.trajectory_points_
);
}
@@ -134,6 +213,8 @@ EDEPTrajectory& EDEPTrajectory::operator=(EDEPTrajectory&& trj) {
std::swap(this->hit_map_, trj.hit_map_);
std::swap(this->entering_map_, trj.entering_map_);
std::swap(this->exiting_map_, trj.exiting_map_);
+ std::swap(this->first_points_, trj.first_points_);
+ std::swap(this->last_points_, trj.last_points_);
std::swap(this->trajectory_points_, trj.trajectory_points_);
for (auto& t:this->children_trajectories_) {
t.parent_trajectory_ = this;
@@ -259,3 +340,18 @@ bool EDEPTrajectory::IsEntering(component component_name) const {
bool EDEPTrajectory::IsExiting(component component_name) const {
return exiting_map_.at(component_name);
}
+
+EDEPTrajectoryHits& EDEPTrajectory::GetFirstPointInDetector(component component_name) {
+
+ if(HasHitInDetector(component_name) && (first_points_.find(component_name) == first_points_.end())){
+ first_points_[component_name].push_back(trajectory_points_.front());
+ }
+ return first_points_.at(component_name);
+}
+
+EDEPTrajectoryHits& EDEPTrajectory::GetLastPointInDetector(component component_name) {
+ if(HasHitInDetector(component_name) && (last_points_.find(component_name) == last_points_.end())){
+ last_points_[component_name].push_back(trajectory_points_.back());
+ }
+ return last_points_.at(component_name);
+}
\ No newline at end of file
--
GitLab
From 830dcae94fae0f889ef24be4e61d670f9e11e9ed Mon Sep 17 00:00:00 2001
From: Riccardo DAmico <rdamico@neutrino-01.novalocal>
Date: Mon, 19 Feb 2024 17:37:12 +0100
Subject: [PATCH 02/12] Corrected if statement in constructor
---
src/EDEPTrajectory.cpp | 8 +++-----
1 file changed, 3 insertions(+), 5 deletions(-)
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index 5034439..9791538 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -99,12 +99,10 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
first_points_[component::GRAIN].push_back(*next_it);
} else if(in_ecal && next_stt){
last_points_[component::ECAL].push_back(*it);
- first_points_[component::STRAW].push_back(*next_it);
-
+ first_points_[component::STRAW].push_back(*next_it);
} else if(in_ecal && next_mag){
last_points_[component::ECAL].push_back(*it);
- first_points_[component::MAGNET].push_back(*next_it);
-
+ first_points_[component::MAGNET].push_back(*next_it);
} else if(in_ecal && next_world){
last_points_[component::ECAL].push_back(*it);
first_points_[component::WORLD].push_back(*next_it);
@@ -133,7 +131,7 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
} else if(in_world && next_ecal){
last_points_[component::WORLD].push_back(*it);
first_points_[component::ECAL].push_back(*next_it);
- } else if(in_world && next_world){
+ } else if(in_world && next_mag){
last_points_[component::WORLD].push_back(*it);
first_points_[component::MAGNET].push_back(*next_it);
}
--
GitLab
From a31e88e9a076adf894c808d26aa1efeef959d4a8 Mon Sep 17 00:00:00 2001
From: Riccardo DAmico <rdamico@neutrino-01.novalocal>
Date: Tue, 24 Sep 2024 17:46:30 +0200
Subject: [PATCH 03/12] Changed sand_inner volume flag from ecal to world names
in EDEPHit. Added generalizations to EDEPTrajectory first and last points.
---
include/EDEPTrajectory.h | 1 -
src/EDEPHit.cpp | 4 +-
src/EDEPTrajectory.cpp | 82 +++++++++++++++++++++++++++++-----------
3 files changed, 61 insertions(+), 26 deletions(-)
diff --git a/include/EDEPTrajectory.h b/include/EDEPTrajectory.h
index 7d8d10c..ae5e4e5 100644
--- a/include/EDEPTrajectory.h
+++ b/include/EDEPTrajectory.h
@@ -76,7 +76,6 @@ class EDEPTrajectory {
bool HasHitBeforeTime(double start_time) const;
bool HasHitAfterTime(double stop_time) const;
bool IsTrajectorySaturated() const;
-
bool Match(std::string volume, std::initializer_list<std::string> names) const;
bool IsEntering(component component_name) const;
bool IsExiting(component component_name) const;
diff --git a/src/EDEPHit.cpp b/src/EDEPHit.cpp
index a3ba6f8..e2073d5 100644
--- a/src/EDEPHit.cpp
+++ b/src/EDEPHit.cpp
@@ -18,6 +18,6 @@ std::map<std::string, component> string_to_component = {
std::initializer_list<std::string> grain_names = {"GRAIN", "GRIAN"};
std::initializer_list<std::string> stt_names = {"horizontalST", "STT", "TrkMod", "CMod", "C3H6Mod"};
-std::initializer_list<std::string> ecal_names = {"ECAL", "kloe_calo_volume", "sand_inner"};
+std::initializer_list<std::string> ecal_names = {"ECAL", "kloe_calo_volume"};
std::initializer_list<std::string> magnet_names = {"KLOE", "Yoke", "Mag"};
-std::initializer_list<std::string> world_names = {"World", "rock", "volSAND"};
+std::initializer_list<std::string> world_names = {"World", "rock", "volSAND" , "volDetEnclosure" , "sand_inner"};
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index 9791538..6a66a1a 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -31,28 +31,55 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
trajectory_points_.push_back(*it);
auto next_it = std::next(it);
-
- if (next_it == trajectory.Points.end()) {
- continue;
- }
-
+
+
auto current_node = GetNode(*it, geo);
- auto current_volume = current_node->GetName();
- auto next_node = GetNode(*next_it, geo);
- auto next_volume = next_node->GetName();
-
+ std::string current_volume;
+ if(current_node != nullptr ){
+ current_volume = current_node->GetName();
+ } else {
+ current_volume = "";
+ }
bool in_grain = Match(current_volume, grain_names);
bool in_stt = Match(current_volume, stt_names);
bool in_ecal = Match(current_volume, ecal_names);
bool in_mag = Match(current_volume, magnet_names);
bool in_world = Match(current_volume, world_names);
+
+ if(it == trajectory.Points.begin()){
+ if(in_ecal) first_points_[component::ECAL].push_back(*it);
+ if(in_grain) first_points_[component::GRAIN].push_back(*it);
+ if(in_stt) first_points_[component::STRAW].push_back(*it);
+ if(in_mag) first_points_[component::MAGNET].push_back(*it);
+ if(in_world) first_points_[component::WORLD].push_back(*it);
+ }
+ if(next_it == trajectory.Points.end()){
+ if(in_ecal) last_points_[component::ECAL].push_back(*it);
+ if(in_grain) last_points_[component::GRAIN].push_back(*it);
+ if(in_stt) last_points_[component::STRAW].push_back(*it);
+ if(in_mag) last_points_[component::MAGNET].push_back(*it);
+ if(in_world) last_points_[component::WORLD].push_back(*it);
+ continue;
+
+ }
+ auto next_node = GetNode(*next_it, geo);
+ std::string next_volume;
+ if(next_node != nullptr){
+ next_volume = next_node->GetName();
+
+ } else {
+ next_volume = "";
+
+ }
+
+
bool next_grain = Match(next_volume, grain_names);
bool next_stt = Match(next_volume, stt_names);
bool next_ecal = Match(next_volume, ecal_names);
bool next_mag = Match(next_volume, magnet_names);
bool next_world = Match(next_volume, world_names);
-
+
if (!exiting_map_[component::GRAIN]) exiting_map_[component::GRAIN] = (in_grain && ( next_stt || next_ecal || next_mag || next_world ));
if (!exiting_map_[component::STRAW]) exiting_map_[component::STRAW] = (in_stt && ( next_grain || next_ecal || next_mag || next_world ));
if (!exiting_map_[component::ECAL]) exiting_map_[component::ECAL] = (in_ecal && ( next_stt || next_grain || next_mag || next_world ));
@@ -64,8 +91,7 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
if (!entering_map_[component::ECAL]) entering_map_[component::ECAL] = (next_ecal && ( in_stt || in_grain || in_mag || in_world ));
if (!entering_map_[component::MAGNET]) entering_map_[component::MAGNET] = (next_mag && ( in_stt || in_ecal || in_grain || in_world ));
if (!entering_map_[component::WORLD]) entering_map_[component::WORLD] = (next_world && ( in_stt || in_ecal || in_mag || in_grain ));
-
-
+
if(in_grain && next_stt){
last_points_[component::GRAIN].push_back(*it);
first_points_[component::STRAW].push_back(*next_it);
@@ -79,7 +105,7 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
last_points_[component::GRAIN].push_back(*it);
first_points_[component::WORLD].push_back(*next_it);
}
-
+
if(in_stt && next_grain){
last_points_[component::STRAW].push_back(*it);
first_points_[component::GRAIN].push_back(*next_it);
@@ -93,7 +119,7 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
last_points_[component::STRAW].push_back(*it);
first_points_[component::WORLD].push_back(*next_it);
}
-
+
if(in_ecal && next_grain){
last_points_[component::ECAL].push_back(*it);
first_points_[component::GRAIN].push_back(*next_it);
@@ -107,7 +133,7 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
last_points_[component::ECAL].push_back(*it);
first_points_[component::WORLD].push_back(*next_it);
}
-
+
if(in_mag && next_grain){
last_points_[component::MAGNET].push_back(*it);
first_points_[component::GRAIN].push_back(*next_it);
@@ -121,7 +147,7 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
last_points_[component::MAGNET].push_back(*it);
first_points_[component::WORLD].push_back(*next_it);
}
-
+
if(in_world && next_grain){
last_points_[component::WORLD].push_back(*it);
first_points_[component::GRAIN].push_back(*next_it);
@@ -134,11 +160,11 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
} else if(in_world && next_mag){
last_points_[component::WORLD].push_back(*it);
first_points_[component::MAGNET].push_back(*next_it);
- }
+ }
}
for (const auto& hmap:hits) {
- std::vector<EDEPHit> edep_hits;
+ std::vector<EDEPHit> edep_hits;
for (uint i = 0; i < hmap.second.size(); i++) {
auto h = hmap.second[i];
if (id_ == h.Contrib[0]) {
@@ -153,6 +179,7 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
}
+
}
EDEPTrajectory::EDEPTrajectory(const EDEPTrajectory& trj) {
@@ -341,15 +368,24 @@ bool EDEPTrajectory::IsExiting(component component_name) const {
EDEPTrajectoryHits& EDEPTrajectory::GetFirstPointInDetector(component component_name) {
- if(HasHitInDetector(component_name) && (first_points_.find(component_name) == first_points_.end())){
- first_points_[component_name].push_back(trajectory_points_.front());
- }
+ if( (IsEntering(component_name)|| HasHitInDetector(component_name)) && (first_points_.find(component_name) == first_points_.end())){
+ first_points_[component_name].push_back(trajectory_points_.front());
+ }// else {
+ // for(int i = 0 ; i< first_points_[component_name].size() ; i ++) {
+ // for(int j = 0 ; j< last_points_[component_name].size() ; j ++){
+ // if(first_points_[component_name].at(i).GetPosition().T() >last_points_[component_name].at(j).GetPosition().T() && i == j){
+ // std::cout << first_points_[component_name].at(i).GetPosition().T() << std::endl;
+ // std::cout << last_points_[component_name].at(j).GetPosition().T() << std::endl;
+ // }
+ // }
+ // }
+ // }
return first_points_.at(component_name);
}
EDEPTrajectoryHits& EDEPTrajectory::GetLastPointInDetector(component component_name) {
- if(HasHitInDetector(component_name) && (last_points_.find(component_name) == last_points_.end())){
- last_points_[component_name].push_back(trajectory_points_.back());
+ if((IsEntering(component_name)|| HasHitInDetector(component_name))&& (last_points_.find(component_name) == last_points_.end())){
+ last_points_[component_name].push_back(trajectory_points_.back());
}
return last_points_.at(component_name);
}
\ No newline at end of file
--
GitLab
From 42ca170ae8278bc78f8851d454d626a3a66af84b Mon Sep 17 00:00:00 2001
From: Valerio Pia <valerio.pia@bo.infn.it>
Date: Sun, 13 Oct 2024 19:52:07 +0200
Subject: [PATCH 04/12] Updated entering/exitig/first/last checks
---
examples/EDEPReader.cpp | 69 +++++-----
include/EDEPTrajectory.h | 9 +-
include/EDEPUtils.h | 3 +
src/EDEPHit.cpp | 5 +
src/EDEPTrajectory.cpp | 284 ++++++++++++++++-----------------------
5 files changed, 169 insertions(+), 201 deletions(-)
diff --git a/examples/EDEPReader.cpp b/examples/EDEPReader.cpp
index c3a6415..01af092 100644
--- a/examples/EDEPReader.cpp
+++ b/examples/EDEPReader.cpp
@@ -22,8 +22,8 @@ int main(int argc, char *argv[])
// auto fMc = TFile::Open("../../mu_pos2_1000_4sett_2.edep.root");
// auto fMc = TFile::Open("/storage/gpfs_data/neutrino/SAND/GRAIN-PHYSICS-CASE/prod/numu_LBNF_GV2_prod/data/numu_LBNF_GV2_prod_1.genie.edep-sim.root");
// auto fMc = TFile::Open("/storage/gpfs_data/neutrino/users/dcasazza/sand-reco/e_pos2_1000.edep.root");
- // auto fMc = TFile::Open("/storage/gpfs_data/neutrino/users/vpia/STTTrackReco/tools/tracklet/DRIFT3.root");
- auto fMc = TFile::Open("/storage/gpfs_data/neutrino/SAND/PRODUCTIONS/PROD/SAND_opt3_DRIFT1_reverse_volSAND/SAND_opt3_DRIFT1_reverse_volSAND_0/sand-drift-spill-revere-volSAND.0.edep.root");
+ auto fMc = TFile::Open("/storage/gpfs_data/neutrino/users/vpia/STTTrackReco/tools/tracklet/DRIFT3.root");
+ // auto fMc = TFile::Open("/storage/gpfs_data/neutrino/SAND/PRODUCTIONS/PROD/SAND_opt3_DRIFT1_reverse_volSAND/SAND_opt3_DRIFT1_reverse_volSAND_0/sand-drift-spill-revere-volSAND.0.edep.root");
// TFile *f = TFile::Open("/storage/gpfs_data/neutrino/users/dcasazza/sand-reco/e_pos2_1000_branch4.digit.Clusters.root");
// TFile *f = TFile::Open("/storage/gpfs_data/neutrino/users/dcasazza/sand-reco/mu_pos2_1000_branch4.digit.Clusters.root");
@@ -37,8 +37,8 @@ int main(int argc, char *argv[])
auto nev = tMc.GetEntries();
// // auto ev_min = std::atoi(argv[1]);
// // auto ev_max = std::atoi(argv[2]);
- auto ev_min = 2;
- auto ev_max = 3;
+ auto ev_min = 1;
+ auto ev_max = 2;
// loop over events
@@ -55,36 +55,45 @@ int main(int argc, char *argv[])
std::string print_string = "";
tree.Print(print_string, 1);
- auto start_point = tree.begin()->GetHitMap().rbegin()->second.rbegin()->GetStart();
- auto stop_point = tree.begin()->GetHitMap().rbegin()->second.rbegin()->GetStop();
- auto middle_point = (start_point + stop_point) * 0.5;
- double distance = 10;
- std::cout << tree.begin()->GetId() << std::endl;
- std::vector<EDEPTrajectory> filteredTrj;
- tree.Filter(std::back_insert_iterator<std::vector<EDEPTrajectory>>(filteredTrj),
- [middle_point, distance](const EDEPTrajectory& trj) { return trj.HasHitNearPoint(middle_point.Vect(), distance);} );
- std::cout << "Number of filtered traj: " << filteredTrj.size() << std::endl;
- for (auto& trj:filteredTrj) {
- std::cout << trj.GetId() << std::endl;
- }
- filteredTrj.clear();
- tree.Filter(std::back_insert_iterator<std::vector<EDEPTrajectory>>(filteredTrj),
- [middle_point, distance](const EDEPTrajectory& trj) { return trj.HasHitNear4DPoint(middle_point, distance, middle_point.T());} );
- std::cout << "Number of filtered traj: " << filteredTrj.size() << std::endl;
- for (auto& trj:filteredTrj) {
- std::cout << trj.GetHitNear4DPoint(middle_point, distance, middle_point.T())->GetId() << std::endl;
- }
+ // auto start_point = tree.begin()->GetHitMap().rbegin()->second.rbegin()->GetStart();
+ // auto stop_point = tree.begin()->GetHitMap().rbegin()->second.rbegin()->GetStop();
+ // auto middle_point = (start_point + stop_point) * 0.5;
+ // double distance = 10;
+ // std::cout << tree.begin()->GetId() << std::endl;
+ // std::vector<EDEPTrajectory> filteredTrj;
+ // tree.Filter(std::back_insert_iterator<std::vector<EDEPTrajectory>>(filteredTrj),
+ // [middle_point, distance](const EDEPTrajectory& trj) { return trj.HasHitNearPoint(middle_point.Vect(), distance);} );
+ // std::cout << "Number of filtered traj: " << filteredTrj.size() << std::endl;
+ // for (auto& trj:filteredTrj) {
+ // std::cout << trj.GetId() << std::endl;
+ // }
+ // filteredTrj.clear();
+ // tree.Filter(std::back_insert_iterator<std::vector<EDEPTrajectory>>(filteredTrj),
+ // [middle_point, distance](const EDEPTrajectory& trj) { return trj.HasHitNear4DPoint(middle_point, distance, middle_point.T());} );
+ // std::cout << "Number of filtered traj: " << filteredTrj.size() << std::endl;
+ // for (auto& trj:filteredTrj) {
+ // std::cout << trj.GetHitNear4DPoint(middle_point, distance, middle_point.T())->GetId() << std::endl;
+ // }
// std::cout << "########### Printing the tree" << std::endl;
// tree.Print(print_string);
- // // LArHit, Straw, EMCalSci
- // component comp = component::ECAL;
- // std::cout << "########### Filtering all the trajectories Entering in " << component_to_string.at(comp) << std::endl;
- // std::vector<EDEPTrajectory> filteredTrj;
- // tree.Filter(std::back_insert_iterator<std::vector<EDEPTrajectory>>(filteredTrj),
- // [comp](const EDEPTrajectory& trj) { return trj.IsEntering(comp);} );
- // std::cout << "Number of filtered traj: " << filteredTrj.size() << std::endl;
+ // LArHit, Straw, EMCalSci, DriftVolume
+ component comp = component::DRIFT;
+ std::cout << "########### Filtering all the trajectories Entering in " << component_to_string.at(comp) << std::endl;
+ std::vector<EDEPTrajectory> filteredTrj;
+ tree.Filter(std::back_insert_iterator<std::vector<EDEPTrajectory>>(filteredTrj),
+ [comp](const EDEPTrajectory& trj) { return trj.IsExiting(comp);} );
+ std::cout << "Number of filtered traj: " << filteredTrj.size() << std::endl;
+
+ for (auto trj:filteredTrj) {
+ for (auto p:trj.GetFirstPointsInDetector(comp)) {
+ p.GetPosition().Print();
+ }
+ for (auto p:trj.GetLastPointsInDetector(comp)) {
+ p.GetPosition().Print();
+ }
+ }
// filteredTrj.clear();
diff --git a/include/EDEPTrajectory.h b/include/EDEPTrajectory.h
index 852c1ac..ec4d97a 100644
--- a/include/EDEPTrajectory.h
+++ b/include/EDEPTrajectory.h
@@ -117,8 +117,8 @@ class EDEPTrajectory {
int GetPDGCode() const {return pdg_code_;};
- EDEPTrajectoryHits& GetFirstPointInDetector(component component_name);
- EDEPTrajectoryHits& GetLastPointInDetector(component component_name);
+ std::vector<EDEPTrajectoryPoint>& GetFirstPointsInDetector(component component_name);
+ std::vector<EDEPTrajectoryPoint>& GetLastPointsInDetector(component component_name);
/**
* @brief Get the initial momentum of this trajectory.
@@ -198,7 +198,6 @@ class EDEPTrajectory {
bool HasHitBeforeTime(double start_time) const;
bool HasHitAfterTime(double stop_time) const;
bool IsTrajectorySaturated() const;
- bool Match(std::string volume, std::initializer_list<std::string> names) const;
bool HasHitInTime(double start_time, double stop_time) const;
bool HasHitWithIdInDetector( int id, component component_name) const;
@@ -237,6 +236,8 @@ class EDEPTrajectory {
bool Match(std::string volume, std::initializer_list<std::string> names) const;
+ void CheckInNext(bool* in, bool* next, TG4TrajectoryPoint it, TG4TrajectoryPoint next_it);
+
friend class EDEPTree;
private:
@@ -246,6 +247,8 @@ class EDEPTrajectory {
std::vector<EDEPTrajectory> children_trajectories_; ///< Children trajectories.
std::map<component, bool> exiting_map_; ///< Map indicating whether the trajectory is exiting a component.
std::map<component, bool> entering_map_; ///< Map indicating whether the trajectory is entering a component.
+ EDEPTrajectoryPoints last_points_; ///< Map of all the first points in each component.
+ EDEPTrajectoryPoints first_points_; ///< Map of all the last points in each component.
EDEPTrajectory* parent_trajectory_; ///< Pointer to the parent trajectory.
int id_; ///< ID of the trajectory.
int parent_id_; ///< Parent ID of the trajectory.
diff --git a/include/EDEPUtils.h b/include/EDEPUtils.h
index 271d094..b323433 100644
--- a/include/EDEPUtils.h
+++ b/include/EDEPUtils.h
@@ -18,6 +18,9 @@ enum class component {
WORLD
};
+extern component components[6];
+
+
extern std::map<component, std::string> component_to_string;
extern std::map<std::string, component> string_to_component;
diff --git a/src/EDEPHit.cpp b/src/EDEPHit.cpp
index 3289519..ada06ce 100644
--- a/src/EDEPHit.cpp
+++ b/src/EDEPHit.cpp
@@ -25,6 +25,11 @@ std::map<std::string, component> string_to_component = {
{ "World", component::WORLD },
};
+/**
+ * @brief List of components in an array.
+ */
+ component components[6] = {component::GRAIN, component::STRAW, component::DRIFT, component::ECAL, component::MAGNET, component::WORLD};
+
/**
* @brief List of string names associated with the GRAIN component.
*/
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index 6c6701a..c54d822 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -22,7 +22,6 @@ TGeoNode* GetNode(const TG4TrajectoryPoint& tpoint, TGeoManager* geo) {
// if (abs(mom.X()) == 0 || abs(mom.Y()) == 0 || abs(mom.Z()) == 0) {
// return node;
// }
- // geo->FindNextBoundary(1000);
geo->FindNextBoundary(1000);
if (geo->GetStep() < 1E-5) {
@@ -33,6 +32,56 @@ TGeoNode* GetNode(const TG4TrajectoryPoint& tpoint, TGeoManager* geo) {
return node;
}
+// Notice: this works (I think) but I don't like it.
+/**
+ * @brief Checks for transitions between detector components in a particle trajectory.
+ *
+ * This function examines the current and next trajectory points to determine
+ * if a transition occurs from one detector component to another.
+ * If a transition is detected, it updates the entering and exiting state
+ * of the respective components and stores the associated trajectory points.
+ *
+ * The `in` and `next` arrays should reflect the following order of components:
+ * - 0: GRAIN
+ * - 1: STRAW
+ * - 2: DRIFT
+ * - 3: ECAL
+ * - 4: MAGNET
+ * - 5: WORLD
+ *
+ * @param in Pointer to a boolean array indicating the components the current
+ * trajectory point is inside.
+ * @param next Pointer to a boolean array indicating the components the next
+ * trajectory point is inside.
+ * @param it The current trajectory point being checked.
+ * @param next_it The next trajectory point being checked.
+ */
+void EDEPTrajectory::CheckInNext(bool* in, bool* next, TG4TrajectoryPoint it, TG4TrajectoryPoint next_it) {
+ for (int i = 0; i < 6; i++) {
+ for (int j = 0; j < 6; j++) {
+ // Notice: this should be fixed in the geometry.
+ if ((components[i] == component::STRAW && components[j] == component::DRIFT) ||
+ (components[j] == component::STRAW && components[i] == component::DRIFT)) {
+ continue;
+ }
+ if (i == j) {
+ continue;
+ }
+ if(in[i] && next[j]) {
+ std::cout << component_to_string[components[i]] << " " << component_to_string[components[j]] << std::endl;
+ last_points_[components[i]].push_back(it);
+ first_points_[components[j]].push_back(next_it);
+
+ exiting_map_[components[i]] = true;
+ entering_map_[components[j]] = true;
+
+ return;
+ }
+ }
+ }
+ return;
+}
+
/**
* @brief Constructs an EDEPTrajectory object from a TG4Trajectory and hit information.
*
@@ -61,163 +110,78 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
break;
}
}
+
+ for (auto comp:string_to_component) {
+ last_points_[comp.second] = {};
+ first_points_[comp.second] = {};
+
+ exiting_map_[comp.second] = false;
+ entering_map_[comp.second] = false;
+ }
for (auto it = trajectory.Points.begin(); it != trajectory.Points.end(); ++it) {
auto next_it = std::next(it);
- auto current_node = GetNode(*it, geo);
- std::string current_volume;
- if(current_node != nullptr ){
- current_volume = current_node->GetName();
- } else {
- current_volume = "";
- }
- bool in_grain = Match(current_volume, grain_names);
- bool in_stt = Match(current_volume, stt_names);
- bool in_ecal = Match(current_volume, ecal_names);
- bool in_mag = Match(current_volume, magnet_names);
- bool in_world = Match(current_volume, world_names);
-
-
- if(it == trajectory.Points.begin()){
- if(in_ecal) first_points_[component::ECAL].push_back(*it);
- if(in_grain) first_points_[component::GRAIN].push_back(*it);
- if(in_stt) first_points_[component::STRAW].push_back(*it);
- if(in_mag) first_points_[component::MAGNET].push_back(*it);
- if(in_world) first_points_[component::WORLD].push_back(*it);
- }
- if(next_it == trajectory.Points.end()){
- if(in_ecal) last_points_[component::ECAL].push_back(*it);
- if(in_grain) last_points_[component::GRAIN].push_back(*it);
- if(in_stt) last_points_[component::STRAW].push_back(*it);
- if(in_mag) last_points_[component::MAGNET].push_back(*it);
- if(in_world) last_points_[component::WORLD].push_back(*it);
- continue;
-
- }
- auto next_node = GetNode(*next_it, geo);
- std::string next_volume;
- if(next_node != nullptr){
- next_volume = next_node->GetName();
- } else {
- next_volume = "";
- }
+ auto current_node = GetNode(*it, geo);
+ std::string current_volume;
+ if(current_node != nullptr ){
+ current_volume = current_node->GetName();
+ } else {
+ current_volume = "";
+ }
- component comp;
- bool in_grain = Match(current_volume, grain_names);
- if(in_grain) comp = component::GRAIN;
-
- bool in_stt = Match(current_volume, stt_names);
- bool in_drift = Match(current_volume, drift_names);
- if(in_stt) comp = component::STRAW;
- if(in_drift) comp = component::DRIFT;
-
- bool in_ecal = Match(current_volume, ecal_names);
- if(in_ecal) comp = component::ECAL;
-
- bool in_mag = Match(current_volume, magnet_names);
- if(in_mag) comp = component::MAGNET;
-
- bool in_world = Match(current_volume, world_names);
- if(in_world) comp = component::WORLD;
-
- // Notice: Ths is useful to debug unincluded volume names
- // if (!(in_grain || in_stt || in_drift || in_ecal || in_mag || in_world)) {
- // std::cout << current_volume << " " << component_to_string[comp] << std::endl;
- // }
- trajectory_points_[comp].push_back(*it);
-
- bool next_grain = Match(next_volume, grain_names);
- bool next_stt = Match(next_volume, stt_names);
- bool next_drift = Match(next_volume, drift_names);
- bool next_ecal = Match(next_volume, ecal_names);
- bool next_mag = Match(next_volume, magnet_names);
- bool next_world = Match(next_volume, world_names);
-
- if (!exiting_map_[component::GRAIN]) exiting_map_[component::GRAIN] = (in_grain && ( next_stt || next_drift || next_ecal || next_mag || next_world ));
- if (!exiting_map_[component::STRAW]) exiting_map_[component::STRAW] = (in_stt && ( next_grain || next_ecal || next_mag || next_world ));
- if (!exiting_map_[component::DRIFT]) exiting_map_[component::DRIFT] = (in_drift && ( next_grain || next_ecal || next_mag || next_world ));
- if (!exiting_map_[component::ECAL]) exiting_map_[component::ECAL] = (in_ecal && ( next_stt || next_drift || next_grain || next_mag || next_world ));
- if (!exiting_map_[component::MAGNET]) exiting_map_[component::MAGNET] = (in_mag && ( next_stt || next_drift || next_ecal || next_grain || next_world ));
- if (!exiting_map_[component::WORLD]) exiting_map_[component::WORLD] = (in_world && ( next_stt || next_drift || next_ecal || next_mag || next_grain ));
-
- if (!entering_map_[component::GRAIN]) entering_map_[component::GRAIN] = (next_grain && ( in_stt || in_drift || in_ecal || in_mag || in_world ));
- if (!entering_map_[component::STRAW]) entering_map_[component::STRAW] = (next_stt && ( in_grain || in_ecal || in_mag || in_world ));
- if (!entering_map_[component::DRIFT]) entering_map_[component::DRIFT] = (next_drift && ( in_grain || in_ecal || in_mag || in_world ));
- if (!entering_map_[component::ECAL]) entering_map_[component::ECAL] = (next_ecal && ( in_stt || in_drift || in_grain || in_mag || in_world ));
- if (!entering_map_[component::MAGNET]) entering_map_[component::MAGNET] = (next_mag && ( in_stt || in_drift || in_ecal || in_grain || in_world ));
- if (!entering_map_[component::WORLD]) entering_map_[component::WORLD] = (next_world && ( in_stt || in_drift || in_ecal || in_mag || in_grain ));
-
-
- if(in_grain && next_stt){
- last_points_[component::GRAIN].push_back(*it);
- first_points_[component::STRAW].push_back(*next_it);
- } else if(in_grain && next_ecal){
- last_points_[component::GRAIN].push_back(*it);
- first_points_[component::ECAL].push_back(*next_it);
- } else if(in_grain && next_mag){
- last_points_[component::GRAIN].push_back(*it);
- first_points_[component::MAGNET].push_back(*next_it);
- } else if(in_grain && next_world){
- last_points_[component::GRAIN].push_back(*it);
- first_points_[component::WORLD].push_back(*next_it);
- }
+ auto next_node = GetNode(*next_it, geo);
+ std::string next_volume;
+ if(next_node != nullptr){
+ next_volume = next_node->GetName();
+ } else {
+ next_volume = "";
+ }
+
+ component comp;
+ bool in_grain = Match(current_volume, grain_names);
+ if(in_grain) comp = component::GRAIN;
- if(in_stt && next_grain){
- last_points_[component::STRAW].push_back(*it);
- first_points_[component::GRAIN].push_back(*next_it);
- } else if(in_stt && next_ecal){
- last_points_[component::STRAW].push_back(*it);
- first_points_[component::ECAL].push_back(*next_it);
- } else if(in_stt && next_mag){
- last_points_[component::STRAW].push_back(*it);
- first_points_[component::MAGNET].push_back(*next_it);
- } else if(in_stt && next_world){
- last_points_[component::STRAW].push_back(*it);
- first_points_[component::WORLD].push_back(*next_it);
- }
+ bool in_stt = Match(current_volume, stt_names);
+ if(in_stt) comp = component::STRAW;
- if(in_ecal && next_grain){
- last_points_[component::ECAL].push_back(*it);
- first_points_[component::GRAIN].push_back(*next_it);
- } else if(in_ecal && next_stt){
- last_points_[component::ECAL].push_back(*it);
- first_points_[component::STRAW].push_back(*next_it);
- } else if(in_ecal && next_mag){
- last_points_[component::ECAL].push_back(*it);
- first_points_[component::MAGNET].push_back(*next_it);
- } else if(in_ecal && next_world){
- last_points_[component::ECAL].push_back(*it);
- first_points_[component::WORLD].push_back(*next_it);
- }
+ bool in_drift = Match(current_volume, drift_names);
+ if(in_drift) comp = component::DRIFT;
- if(in_mag && next_grain){
- last_points_[component::MAGNET].push_back(*it);
- first_points_[component::GRAIN].push_back(*next_it);
- } else if(in_mag && next_stt){
- last_points_[component::MAGNET].push_back(*it);
- first_points_[component::STRAW].push_back(*next_it);
- } else if(in_mag && next_ecal){
- last_points_[component::MAGNET].push_back(*it);
- first_points_[component::ECAL].push_back(*next_it);
- } else if(in_mag && next_world){
- last_points_[component::MAGNET].push_back(*it);
- first_points_[component::WORLD].push_back(*next_it);
- }
+ bool in_ecal = Match(current_volume, ecal_names);
+ if(in_ecal) comp = component::ECAL;
- if(in_world && next_grain){
- last_points_[component::WORLD].push_back(*it);
- first_points_[component::GRAIN].push_back(*next_it);
- } else if(in_world && next_stt){
- last_points_[component::WORLD].push_back(*it);
- first_points_[component::STRAW].push_back(*next_it);
- } else if(in_world && next_ecal){
- last_points_[component::WORLD].push_back(*it);
- first_points_[component::ECAL].push_back(*next_it);
- } else if(in_world && next_mag){
- last_points_[component::WORLD].push_back(*it);
- first_points_[component::MAGNET].push_back(*next_it);
- }
+ bool in_mag = Match(current_volume, magnet_names);
+ if(in_mag) comp = component::MAGNET;
+
+ bool in_world = Match(current_volume, world_names);
+ if(in_world) comp = component::WORLD;
+
+ // Notice: Ths is useful to debug unincluded volume names
+ // if (!(in_grain || in_stt || in_drift || in_ecal || in_mag || in_world)) {
+ // std::cout << current_volume << " " << component_to_string[comp] << std::endl;
+ // }
+ trajectory_points_[comp].push_back(*it);
+
+ if(it == trajectory.Points.begin()){
+ first_points_[comp].push_back(*it);
+ }
+ if(next_it == trajectory.Points.end()){
+ last_points_[comp].push_back(*it);
+ continue;
+ }
+
+ bool next_grain = Match(next_volume, grain_names);
+ bool next_stt = Match(next_volume, stt_names);
+ bool next_drift = Match(next_volume, drift_names);
+ bool next_ecal = Match(next_volume, ecal_names);
+ bool next_mag = Match(next_volume, magnet_names);
+ bool next_world = Match(next_volume, world_names);
+
+ bool in[6] = {in_grain, in_stt, in_drift, in_ecal, in_mag, in_world};
+ bool next[6] = {next_grain, next_stt, next_drift, next_ecal, next_mag, next_world};
+ CheckInNext(in, next, *it, *next_it);
+
}
for (const auto& hmap:hits) {
@@ -268,7 +232,7 @@ EDEPTrajectory& EDEPTrajectory::operator=(const EDEPTrajectory& trj) {
this->hit_map_ = trj.hit_map_;
this->entering_map_ = trj.entering_map_;
this->exiting_map_ = trj.exiting_map_;
- this -> last_points_ = trj.last_points_;
+ this ->last_points_ = trj.last_points_;
this-> first_points_ = trj.first_points_;
this->trajectory_points_ = trj.trajectory_points_;
return *this;
@@ -599,26 +563,10 @@ bool EDEPTrajectory::IsExiting(component component_name) const {
return exiting_map_.at(component_name);
}
-EDEPTrajectoryHits& EDEPTrajectory::GetFirstPointInDetector(component component_name) {
-
- if( (IsEntering(component_name)|| HasHitInDetector(component_name)) && (first_points_.find(component_name) == first_points_.end())){
- first_points_[component_name].push_back(trajectory_points_.front());
- }// else {
- // for(int i = 0 ; i< first_points_[component_name].size() ; i ++) {
- // for(int j = 0 ; j< last_points_[component_name].size() ; j ++){
- // if(first_points_[component_name].at(i).GetPosition().T() >last_points_[component_name].at(j).GetPosition().T() && i == j){
- // std::cout << first_points_[component_name].at(i).GetPosition().T() << std::endl;
- // std::cout << last_points_[component_name].at(j).GetPosition().T() << std::endl;
- // }
- // }
- // }
- // }
+std::vector<EDEPTrajectoryPoint>& EDEPTrajectory::GetFirstPointsInDetector(component component_name) {
return first_points_.at(component_name);
}
-EDEPTrajectoryHits& EDEPTrajectory::GetLastPointInDetector(component component_name) {
- if((IsEntering(component_name)|| HasHitInDetector(component_name))&& (last_points_.find(component_name) == last_points_.end())){
- last_points_[component_name].push_back(trajectory_points_.back());
- }
+std::vector<EDEPTrajectoryPoint>& EDEPTrajectory::GetLastPointsInDetector(component component_name) {
return last_points_.at(component_name);
}
\ No newline at end of file
--
GitLab
From 9fd6bf1b7f91f09ca0e38450e9295bcde65e189a Mon Sep 17 00:00:00 2001
From: Valerio Pia <vpia@bo.infn.it>
Date: Sun, 13 Oct 2024 19:58:20 +0200
Subject: [PATCH 05/12] Simple test
---
examples/EDEPReader.cpp | 32 ++++++++++++++++----------------
1 file changed, 16 insertions(+), 16 deletions(-)
diff --git a/examples/EDEPReader.cpp b/examples/EDEPReader.cpp
index 01af092..b13f178 100644
--- a/examples/EDEPReader.cpp
+++ b/examples/EDEPReader.cpp
@@ -78,22 +78,22 @@ int main(int argc, char *argv[])
// std::cout << "########### Printing the tree" << std::endl;
// tree.Print(print_string);
- // LArHit, Straw, EMCalSci, DriftVolume
- component comp = component::DRIFT;
- std::cout << "########### Filtering all the trajectories Entering in " << component_to_string.at(comp) << std::endl;
- std::vector<EDEPTrajectory> filteredTrj;
- tree.Filter(std::back_insert_iterator<std::vector<EDEPTrajectory>>(filteredTrj),
- [comp](const EDEPTrajectory& trj) { return trj.IsExiting(comp);} );
- std::cout << "Number of filtered traj: " << filteredTrj.size() << std::endl;
-
- for (auto trj:filteredTrj) {
- for (auto p:trj.GetFirstPointsInDetector(comp)) {
- p.GetPosition().Print();
- }
- for (auto p:trj.GetLastPointsInDetector(comp)) {
- p.GetPosition().Print();
- }
- }
+ // // LArHit, Straw, EMCalSci, DriftVolume
+ // component comp = component::DRIFT;
+ // std::cout << "########### Filtering all the trajectories Entering in " << component_to_string.at(comp) << std::endl;
+ // std::vector<EDEPTrajectory> filteredTrj;
+ // tree.Filter(std::back_insert_iterator<std::vector<EDEPTrajectory>>(filteredTrj),
+ // [comp](const EDEPTrajectory& trj) { return trj.IsExiting(comp);} );
+ // std::cout << "Number of filtered traj: " << filteredTrj.size() << std::endl;
+
+ // for (auto trj:filteredTrj) {
+ // for (auto p:trj.GetFirstPointsInDetector(comp)) {
+ // p.GetPosition().Print();
+ // }
+ // for (auto p:trj.GetLastPointsInDetector(comp)) {
+ // p.GetPosition().Print();
+ // }
+ // }
// filteredTrj.clear();
--
GitLab
From ee33ef0525d9d0b5961a6897099d559c4b98d3b0 Mon Sep 17 00:00:00 2001
From: Valerio Pia <vpia@bo.infn.it>
Date: Mon, 14 Oct 2024 10:44:35 +0200
Subject: [PATCH 06/12] Removed cout
---
src/EDEPTrajectory.cpp | 1 -
1 file changed, 1 deletion(-)
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index c54d822..a415066 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -68,7 +68,6 @@ void EDEPTrajectory::CheckInNext(bool* in, bool* next, TG4TrajectoryPoint it, TG
continue;
}
if(in[i] && next[j]) {
- std::cout << component_to_string[components[i]] << " " << component_to_string[components[j]] << std::endl;
last_points_[components[i]].push_back(it);
first_points_[components[j]].push_back(next_it);
--
GitLab
From 59d6652699f6d542303fe69a23ef85ac58f9724f Mon Sep 17 00:00:00 2001
From: Valerio Pia <vpia@bo.infn.it>
Date: Tue, 15 Oct 2024 14:40:55 +0200
Subject: [PATCH 07/12] Added temporary handling of drift/stt geometries
---
src/EDEPTrajectory.cpp | 35 ++++++++++++++++++++++++++++-------
1 file changed, 28 insertions(+), 7 deletions(-)
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index a415066..331c424 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -141,12 +141,6 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
bool in_grain = Match(current_volume, grain_names);
if(in_grain) comp = component::GRAIN;
- bool in_stt = Match(current_volume, stt_names);
- if(in_stt) comp = component::STRAW;
-
- bool in_drift = Match(current_volume, drift_names);
- if(in_drift) comp = component::DRIFT;
-
bool in_ecal = Match(current_volume, ecal_names);
if(in_ecal) comp = component::ECAL;
@@ -156,11 +150,38 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
bool in_world = Match(current_volume, world_names);
if(in_world) comp = component::WORLD;
+
+ // Notice: this should simply be 'trajectory_points_[comp].push_back(*it);'
+ // but stt and drift share some names and only the last one will be added
+ // must be fixed in the geometry
+ if (in_grain || in_ecal || in_mag || in_world) {
+ trajectory_points_[comp].push_back(*it);
+ }
+
+ bool in_stt = Match(current_volume, stt_names);
+ bool in_drift = Match(current_volume, drift_names);
+
+ if(in_stt ) comp = component::STRAW;
+ if(in_drift) comp = component::DRIFT;
+
+
+ if (in_stt && in_drift) {
+ trajectory_points_[component::STRAW].push_back(*it);
+ trajectory_points_[component::DRIFT].push_back(*it);
+ }
+ if (!in_stt && in_drift) {
+ trajectory_points_[component::DRIFT].push_back(*it);
+ }
+ if (in_stt && !in_drift) {
+ trajectory_points_[component::STRAW].push_back(*it);
+ }
+
+
+
// Notice: Ths is useful to debug unincluded volume names
// if (!(in_grain || in_stt || in_drift || in_ecal || in_mag || in_world)) {
// std::cout << current_volume << " " << component_to_string[comp] << std::endl;
// }
- trajectory_points_[comp].push_back(*it);
if(it == trajectory.Points.begin()){
first_points_[comp].push_back(*it);
--
GitLab
From e5a61e18cbaec353a9453d07802a45fa8f147848 Mon Sep 17 00:00:00 2001
From: Valerio Pia <vpia@bo.infn.it>
Date: Sat, 19 Oct 2024 00:57:57 +0200
Subject: [PATCH 08/12] Fixed bug on when to check for next_it
---
examples/EDEPReader.cpp | 8 +++++---
src/EDEPTrajectory.cpp | 41 +++++++++++++++++++----------------------
2 files changed, 24 insertions(+), 25 deletions(-)
diff --git a/examples/EDEPReader.cpp b/examples/EDEPReader.cpp
index b13f178..fddd97b 100644
--- a/examples/EDEPReader.cpp
+++ b/examples/EDEPReader.cpp
@@ -22,7 +22,9 @@ int main(int argc, char *argv[])
// auto fMc = TFile::Open("../../mu_pos2_1000_4sett_2.edep.root");
// auto fMc = TFile::Open("/storage/gpfs_data/neutrino/SAND/GRAIN-PHYSICS-CASE/prod/numu_LBNF_GV2_prod/data/numu_LBNF_GV2_prod_1.genie.edep-sim.root");
// auto fMc = TFile::Open("/storage/gpfs_data/neutrino/users/dcasazza/sand-reco/e_pos2_1000.edep.root");
- auto fMc = TFile::Open("/storage/gpfs_data/neutrino/users/vpia/STTTrackReco/tools/tracklet/DRIFT3.root");
+ // auto fMc = TFile::Open("/storage/gpfs_data/neutrino/users/vpia/STTTrackReco/tools/tracklet/DRIFT3.root");
+ // auto fMc = TFile::Open("/storage/gpfs_data/neutrino/users/rdamico/lotsofnumu/production_14/events-in-rockBox_lv.4.edep.root");
+ auto fMc = TFile::Open("/storage/gpfs_data/neutrino/users/dcasazza/lot-of-nue/production_merge30cm/events-in-STTtracker.33.edep.root");
// auto fMc = TFile::Open("/storage/gpfs_data/neutrino/SAND/PRODUCTIONS/PROD/SAND_opt3_DRIFT1_reverse_volSAND/SAND_opt3_DRIFT1_reverse_volSAND_0/sand-drift-spill-revere-volSAND.0.edep.root");
// TFile *f = TFile::Open("/storage/gpfs_data/neutrino/users/dcasazza/sand-reco/e_pos2_1000_branch4.digit.Clusters.root");
// TFile *f = TFile::Open("/storage/gpfs_data/neutrino/users/dcasazza/sand-reco/mu_pos2_1000_branch4.digit.Clusters.root");
@@ -37,8 +39,8 @@ int main(int argc, char *argv[])
auto nev = tMc.GetEntries();
// // auto ev_min = std::atoi(argv[1]);
// // auto ev_max = std::atoi(argv[2]);
- auto ev_min = 1;
- auto ev_max = 2;
+ auto ev_min = 0;
+ auto ev_max = 250;
// loop over events
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index 331c424..24609fb 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -129,14 +129,6 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
current_volume = "";
}
- auto next_node = GetNode(*next_it, geo);
- std::string next_volume;
- if(next_node != nullptr){
- next_volume = next_node->GetName();
- } else {
- next_volume = "";
- }
-
component comp;
bool in_grain = Match(current_volume, grain_names);
if(in_grain) comp = component::GRAIN;
@@ -150,6 +142,17 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
bool in_world = Match(current_volume, world_names);
if(in_world) comp = component::WORLD;
+ bool in_stt = Match(current_volume, stt_names);
+ bool in_drift = Match(current_volume, drift_names);
+
+ if(in_stt ) comp = component::STRAW;
+ if(in_drift) comp = component::DRIFT;
+
+ // Notice: Ths is useful to debug unincluded volume names
+ // if (!(in_grain || in_stt || in_drift || in_ecal || in_mag || in_world)) {
+ // continue;
+ // std::cout << current_volume << " " << component_to_string[comp] << std::endl;
+ // }
// Notice: this should simply be 'trajectory_points_[comp].push_back(*it);'
// but stt and drift share some names and only the last one will be added
@@ -158,13 +161,6 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
trajectory_points_[comp].push_back(*it);
}
- bool in_stt = Match(current_volume, stt_names);
- bool in_drift = Match(current_volume, drift_names);
-
- if(in_stt ) comp = component::STRAW;
- if(in_drift) comp = component::DRIFT;
-
-
if (in_stt && in_drift) {
trajectory_points_[component::STRAW].push_back(*it);
trajectory_points_[component::DRIFT].push_back(*it);
@@ -176,13 +172,6 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
trajectory_points_[component::STRAW].push_back(*it);
}
-
-
- // Notice: Ths is useful to debug unincluded volume names
- // if (!(in_grain || in_stt || in_drift || in_ecal || in_mag || in_world)) {
- // std::cout << current_volume << " " << component_to_string[comp] << std::endl;
- // }
-
if(it == trajectory.Points.begin()){
first_points_[comp].push_back(*it);
}
@@ -191,6 +180,14 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
continue;
}
+ auto next_node = GetNode(*next_it, geo);
+ std::string next_volume;
+ if(next_node != nullptr){
+ next_volume = next_node->GetName();
+ } else {
+ next_volume = "";
+ }
+
bool next_grain = Match(next_volume, grain_names);
bool next_stt = Match(next_volume, stt_names);
bool next_drift = Match(next_volume, drift_names);
--
GitLab
From 204779ca699b211455f6718fc830ed3db9279752 Mon Sep 17 00:00:00 2001
From: Valerio Pia <vpia@bo.infn.it>
Date: Sat, 19 Oct 2024 00:59:22 +0200
Subject: [PATCH 09/12] Moved hit_map creation in CreateTree to avoid looping
all hits every trj
---
include/EDEPTrajectory.h | 2 +-
src/EDEPTrajectory.cpp | 20 +++++---------------
src/EDEPTree.cpp | 9 ++++++++-
3 files changed, 14 insertions(+), 17 deletions(-)
diff --git a/include/EDEPTrajectory.h b/include/EDEPTrajectory.h
index ec4d97a..f440898 100644
--- a/include/EDEPTrajectory.h
+++ b/include/EDEPTrajectory.h
@@ -29,7 +29,7 @@ class EDEPTrajectory {
parent_id_(trajectory.GetParentId()),
pdg_code_(trajectory.GetPDGCode()) {};
- EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegmentDetectors& hits, const TG4PrimaryVertexContainer& primaries, TGeoManager* geo);
+ EDEPTrajectory(const TG4Trajectory& trajectory, const std::map<int, std::map<component, std::vector<EDEPHit>>>& hit_map, const TG4PrimaryVertexContainer& primaries, TGeoManager* geo);
EDEPTrajectory(const EDEPTrajectory& trj);
EDEPTrajectory(EDEPTrajectory&& trj);
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index 24609fb..64297b9 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -94,7 +94,7 @@ void EDEPTrajectory::CheckInNext(bool* in, bool* next, TG4TrajectoryPoint it, TG
* @param primaries A container for primary vertex information to track interaction number and reaction type.
* @param geo The TGeoManager object representing the detector geometry for identifying volume regions.
*/
-EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegmentDetectors& hits, const TG4PrimaryVertexContainer& primaries, TGeoManager* geo) :
+EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const std::map<int, std::map<component, std::vector<EDEPHit>>>& hit_map, const TG4PrimaryVertexContainer& primaries, TGeoManager* geo) :
p0_(trajectory.GetInitialMomentum()),
parent_trajectory_(nullptr),
id_(trajectory.GetTrackId()),
@@ -201,23 +201,13 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
}
- for (const auto& hmap:hits) {
- std::vector<EDEPHit> edep_hits;
- for (uint i = 0; i < hmap.second.size(); i++) {
- auto h = hmap.second[i];
- if (id_ == h.Contrib[0]) {
- edep_hits.push_back(EDEPHit(h, i));
- }
- }
- if (!edep_hits.empty()) {
- std::sort(edep_hits.begin(), edep_hits.end(),
+ if (hit_map.find(id_) != hit_map.end()) {
+ hit_map_ = hit_map.at(id_);
+ for (auto& hits:hit_map_) {
+ std::sort(hits.second.begin(), hits.second.end(),
[](EDEPHit i,EDEPHit j) { return i.GetStart().T() < j.GetStart().T();});
- hit_map_[string_to_component[hmap.first]] = edep_hits;
}
}
-
-
-
}
/**
diff --git a/src/EDEPTree.cpp b/src/EDEPTree.cpp
index 07496c6..544b48d 100644
--- a/src/EDEPTree.cpp
+++ b/src/EDEPTree.cpp
@@ -126,9 +126,16 @@ void EDEPTree::CreateTree(const std::vector<EDEPTrajectory>& trajectories_vect)
*/
void EDEPTree::InizializeFromEdep(const TG4Event& edep_event, TGeoManager* geo) {
this->GetChildrenTrajectories().clear();
+ std::map<int, std::map<component, std::vector<EDEPHit>>> hit_map;
+ for (auto hmap:edep_event.SegmentDetectors) {
+ for (uint i = 0; i < hmap.second.size(); i++) {
+ auto h = hmap.second[i];
+ hit_map[h.Contrib[0]][string_to_component[hmap.first]].push_back(EDEPHit(h, i));
+ }
+ }
std::vector<EDEPTrajectory> trajectories;
for (auto trj:edep_event.Trajectories) {
- trajectories.push_back(EDEPTrajectory(trj, edep_event.SegmentDetectors, edep_event.Primaries, geo));
+ trajectories.push_back(EDEPTrajectory(trj, hit_map, edep_event.Primaries, geo));
}
CreateTree(trajectories);
}
--
GitLab
From 7288d80619ec5f961b1d7aa8bf86468714ee0092 Mon Sep 17 00:00:00 2001
From: Valerio Pia <vpia@bo.infn.it>
Date: Wed, 23 Oct 2024 10:45:57 +0200
Subject: [PATCH 10/12] Added getters of all trj point in a sorted vector
---
include/EDEPTrajectory.h | 6 ++++++
src/EDEPTrajectory.cpp | 12 ++++++++++++
2 files changed, 18 insertions(+)
diff --git a/include/EDEPTrajectory.h b/include/EDEPTrajectory.h
index f440898..3cc2deb 100644
--- a/include/EDEPTrajectory.h
+++ b/include/EDEPTrajectory.h
@@ -150,6 +150,12 @@ class EDEPTrajectory {
*/
const EDEPTrajectoryPoints& GetTrajectoryPoints() const {return trajectory_points_;};
+ /**
+ * @brief Get the trajectory points associated with this trajectory.
+ * @return All trajectory points associated with this trajectory ordered by increasing times.
+ */
+ std::vector<EDEPTrajectoryPoint> GetTrajectoryPointsVect();
+
// Setters
/**
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index 64297b9..aa6540d 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -576,4 +576,16 @@ std::vector<EDEPTrajectoryPoint>& EDEPTrajectory::GetFirstPointsInDetector(compo
std::vector<EDEPTrajectoryPoint>& EDEPTrajectory::GetLastPointsInDetector(component component_name) {
return last_points_.at(component_name);
+}
+
+std::vector<EDEPTrajectoryPoint> EDEPTrajectory::GetTrajectoryPointsVect() {
+ std::vector<EDEPTrajectoryPoint> points;
+ for (auto& comp:trajectory_points_) {
+ points.insert(points.end(), comp.second.begin(), comp.second.end());
+ }
+
+ std::sort(points.begin(), points.end(),
+ [](EDEPTrajectoryPoint i,EDEPTrajectoryPoint j) { return i.GetPosition().T() < j.GetPosition().T();});
+
+ return points;
}
\ No newline at end of file
--
GitLab
From 0750173dc614bc084d9c7dcc3bb222815baccfec Mon Sep 17 00:00:00 2001
From: dcasazza <denise.casazza@fe.infn.it>
Date: Thu, 24 Oct 2024 16:32:09 +0200
Subject: [PATCH 11/12] Added operator== for trajectoryPoint
---
include/EDEPTrajectoryPoint.h | 9 +++++++++
1 file changed, 9 insertions(+)
diff --git a/include/EDEPTrajectoryPoint.h b/include/EDEPTrajectoryPoint.h
index d524d2c..f73a73f 100644
--- a/include/EDEPTrajectoryPoint.h
+++ b/include/EDEPTrajectoryPoint.h
@@ -24,6 +24,15 @@ class EDEPTrajectoryPoint {
*/
~EDEPTrajectoryPoint() {};
+ bool operator==(const EDEPTrajectoryPoint& p) {
+ return (
+ this->position_ == p.GetPosition() &&
+ this->momentum_ == p.GetMomentum() &&
+ this->process_ == p.GetProcess() &&
+ this->sub_process_ == p.GetSubprocess()
+ );
+ };
+
/**
* @brief Get the position of the trajectory point.
* @return The position as a TLorentzVector.
--
GitLab
From 9de1d346ea609eff63ba5c0f9f318ece2e085cc0 Mon Sep 17 00:00:00 2001
From: Valerio Pia <vpia@bo.infn.it>
Date: Sun, 27 Oct 2024 17:58:06 +0100
Subject: [PATCH 12/12] Fix for particles at rest
---
src/EDEPTrajectory.cpp | 11 ++++++-----
1 file changed, 6 insertions(+), 5 deletions(-)
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index aa6540d..abc129d 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -17,11 +17,12 @@ TGeoNode* GetNode(const TG4TrajectoryPoint& tpoint, TGeoManager* geo) {
geo->SetCurrentDirection(mom.X(), mom.Y(), mom.Z());
- // Notice: this was added for the KF. Check if really needed.
- // Commented at the moment.
- // if (abs(mom.X()) == 0 || abs(mom.Y()) == 0 || abs(mom.Z()) == 0) {
- // return node;
- // }
+ // Notice: this is a fix for cases of P = 0.
+ // Particle not moving means FindNextBoundary will not converge
+ // Is this edepsim?
+ if (abs(mom.X()) == 0 && abs(mom.Y()) == 0 && abs(mom.Z()) == 0) {
+ return node;
+ }
geo->FindNextBoundary(1000);
if (geo->GetStep() < 1E-5) {
--
GitLab