diff --git a/include/EDEPTrajectory.h b/include/EDEPTrajectory.h
index 87e4a44409e310fab598734b14adb4f50d484e1c..852c1ac1630cd036a5573b866e15ebbdda4fabe1 100644
--- a/include/EDEPTrajectory.h
+++ b/include/EDEPTrajectory.h
@@ -115,6 +115,10 @@ class EDEPTrajectory {
* @return The PDG code of this trajectory.
*/
int GetPDGCode() const {return pdg_code_;};
+
+
+ EDEPTrajectoryHits& GetFirstPointInDetector(component component_name);
+ EDEPTrajectoryHits& GetLastPointInDetector(component component_name);
/**
* @brief Get the initial momentum of this trajectory.
@@ -194,6 +198,7 @@ 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;
diff --git a/src/EDEPHit.cpp b/src/EDEPHit.cpp
index 3c79c79cf63a85f7897a68e2856407b4a5bdb468..3289519868386af7b76887b6db50a8f5ab6fd887 100644
--- a/src/EDEPHit.cpp
+++ b/src/EDEPHit.cpp
@@ -53,5 +53,4 @@ std::initializer_list<std::string> magnet_names = {"KLOE", "Yoke", "Mag"};
/**
* @brief List of string names associated with the WORLD component.
*/
-std::initializer_list<std::string> world_names = {"World", "rock", "volSAND", "sand_inner", "Enclosure"};
-
+std::initializer_list<std::string> world_names = {"World", "rock", "volSAND", "sand_inner", "Enclosure", "volDetEnclosure"};
diff --git a/src/EDEPTrajectory.cpp b/src/EDEPTrajectory.cpp
index da6db3a1ecce9a9b0cb3f9742aa290acaea3bd26..6c6701a1e9cfadbe23c50d03deef12486cfcc68f 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -65,63 +65,163 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
for (auto it = trajectory.Points.begin(); it != trajectory.Points.end(); ++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();
-
- component comp;
- bool in_grain = Match(current_volume, grain_names);
- if(in_grain) comp = component::GRAIN;
+ 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 = "";
+ }
- 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;
+ 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);
+ }
- bool in_ecal = Match(current_volume, ecal_names);
- if(in_ecal) comp = component::ECAL;
+ 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_mag = Match(current_volume, magnet_names);
- if(in_mag) comp = component::MAGNET;
+ 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_world = Match(current_volume, world_names);
- if(in_world) comp = component::WORLD;
+ 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);
+ }
- // 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_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);
+ }
}
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]) {
@@ -134,6 +234,9 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
hit_map_[string_to_component[hmap.first]] = edep_hits;
}
}
+
+
+
}
/**
@@ -165,7 +268,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;
}
@@ -182,6 +287,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_
);
}
@@ -213,6 +320,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;
@@ -489,3 +598,27 @@ 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( (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((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