diff --git a/examples/EDEPReader.cpp b/examples/EDEPReader.cpp
index c3a6415b713b52b8c2ebca3b6a14898aba318889..fddd97b857bbe57918994a5cd7ea1d2fe31e92ea 100644
--- a/examples/EDEPReader.cpp
+++ b/examples/EDEPReader.cpp
@@ -23,7 +23,9 @@ int main(int argc, char *argv[])
// 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/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 = 2;
- auto ev_max = 3;
+ auto ev_min = 0;
+ auto ev_max = 250;
// loop over events
@@ -55,36 +57,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;
+ // // 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.IsEntering(comp);} );
+ // [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 87e4a44409e310fab598734b14adb4f50d484e1c..3cc2deb6073813ed8c54eed191646f700d5201a2 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);
@@ -115,6 +115,10 @@ class EDEPTrajectory {
* @return The PDG code of this trajectory.
*/
int GetPDGCode() const {return pdg_code_;};
+
+
+ std::vector<EDEPTrajectoryPoint>& GetFirstPointsInDetector(component component_name);
+ std::vector<EDEPTrajectoryPoint>& GetLastPointsInDetector(component component_name);
/**
* @brief Get the initial momentum of this trajectory.
@@ -146,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
/**
@@ -232,6 +242,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:
@@ -241,6 +253,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/EDEPTrajectoryPoint.h b/include/EDEPTrajectoryPoint.h
index d524d2c3ea3e7f36e595cf23ead639bb36161308..f73a73f33413131a3c9c6e08f65a97e47b109b16 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.
diff --git a/include/EDEPUtils.h b/include/EDEPUtils.h
index 271d094d29b356080bca08f5df7aad03e0a6dd78..b323433ca855e107a1c9cfa07cf178ceeb2e0b30 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 3c79c79cf63a85f7897a68e2856407b4a5bdb468..ada06ce2aca7dfa7e38fdec7f9f4616ac59ea449 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.
*/
@@ -53,5 +58,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..abc129d7ab0d42f7a14e92e9026ea4f77bb0dbf0 100644
--- a/src/EDEPTrajectory.cpp
+++ b/src/EDEPTrajectory.cpp
@@ -17,12 +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;
- // }
- // geo->FindNextBoundary(1000);
+ // 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) {
@@ -33,6 +33,55 @@ 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]) {
+ 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.
*
@@ -46,7 +95,7 @@ TGeoNode* GetNode(const TG4TrajectoryPoint& tpoint, TGeoManager* geo) {
* @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()),
@@ -61,28 +110,30 @@ 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);
- 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 = "";
+ }
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;
@@ -92,12 +143,52 @@ 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)) {
- // std::cout << current_volume << " " << component_to_string[comp] << std::endl;
+ // continue;
+ // std::cout << current_volume << " " << component_to_string[comp] << std::endl;
// }
- trajectory_points_[comp].push_back(*it);
-
+
+ // 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);
+ }
+
+ 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);
+ }
+
+ if(it == trajectory.Points.begin()){
+ first_points_[comp].push_back(*it);
+ }
+ if(next_it == trajectory.Points.end()){
+ last_points_[comp].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_drift = Match(next_volume, drift_names);
@@ -105,33 +196,17 @@ EDEPTrajectory::EDEPTrajectory(const TG4Trajectory& trajectory, const TG4HitSegm
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 ));
+ 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) {
- 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;
}
}
}
@@ -165,7 +240,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 +259,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 +292,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 +570,23 @@ bool EDEPTrajectory::IsEntering(component component_name) const {
bool EDEPTrajectory::IsExiting(component component_name) const {
return exiting_map_.at(component_name);
}
+
+std::vector<EDEPTrajectoryPoint>& EDEPTrajectory::GetFirstPointsInDetector(component component_name) {
+ return first_points_.at(component_name);
+}
+
+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
diff --git a/src/EDEPTree.cpp b/src/EDEPTree.cpp
index 07496c6c95c03001cf676cab6eaa2337b704a00e..544b48d7ce3b6c8df25e74201279cfe4b7351960 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);
}