Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix true visible energy #152

Merged
merged 2 commits into from
Sep 18, 2024
Merged

Fix true visible energy #152

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
345b234
Do not add reco hits that are not in the TMS. Fixes true visible ener…
jdkio Sep 11, 2024
f7dc056
Fix warning to only say something if the map is non-zero size
jdkio Sep 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 25 additions & 43 deletions src/TMS_Event.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,9 @@ TMS_Event::TMS_Event() {


static bool TMS_TrueParticle_NotWorthSaving(TMS_TrueParticle tp) {
//if (tp.GetTrueVisibleEnergy() == 0 && !tp.IsPrimary()) return true;
if (!tp.IsPrimary()) return true;
if (tp.GetTrueVisibleEnergy() == 0 && !tp.IsPrimary()) return true;
// Don't worry about really low visible energy
if (tp.GetTrueVisibleEnergy() < 0.5 && !tp.IsPrimary()) return true;
else return false;
};

Expand All @@ -33,7 +34,7 @@ void TMS_Event::ProcessTG4Event(TG4Event &event, bool FillEvent) {
bool TMSOnly = false;
bool TMSLArOnly = false;
bool OnlyPrimary = false;
bool OnlyPrimaryOrInteresting = true;
bool OnlyPrimaryOrInteresting = false;
bool LightWeight = TMS_Manager::GetInstance().Get_LightWeight_Truth();

int nPrimary = 0;
Expand Down Expand Up @@ -252,26 +253,18 @@ void TMS_Event::ProcessTG4Event(TG4Event &event, bool FillEvent) {
}
else vertex_id = value->second;
TMS_Hit hit = TMS_Hit(edep_hit, vertex_id);
TMS_Hits.push_back(std::move(hit));

// todo, maybe skip for michel electrons or late neutrons
for (size_t i = 0; i < hit.GetTrueHit().GetNTrueParticles(); i++) {
TrueVisibleEnergyPerVertex[hit.GetTrueHit().GetVertexIds(i)] += hit.GetTrueHit().GetEnergyShare(i);
TrueVisibleEnergyPerParticle[hit.GetTrueHit().GetPrimaryIds(i)] += hit.GetTrueHit().GetEnergyShare(i);
}

// Loop through the True Particle list and associate
/*
// Now associate the hits with the muon
int PrimaryId = edep_hit.GetPrimaryId();
for (auto &TrueParticle : TMS_TrueParticles) {
// Check the primary ID
if (TrueParticle.GetTrackId() != PrimaryId) continue;
TLorentzVector Position = (edep_hit.GetStop()+edep_hit.GetStart());
Position *= 0.5;
TrueParticle.AddPoint(Position);
int barnum = hit.GetBarNumber();
// Only add if within the TMS
// Can't use x,y or z because geometry might change. But we know things aren't set if there's no bar number
if (barnum >= 0) {
TMS_Hits.push_back(std::move(hit));

// todo, maybe skip for michel electrons or late neutrons
for (size_t i = 0; i < hit.GetTrueHit().GetNTrueParticles(); i++) {
TrueVisibleEnergyPerVertex[hit.GetTrueHit().GetVertexIds(i)] += hit.GetTrueHit().GetEnergyShare(i);
TrueVisibleEnergyPerParticle[hit.GetTrueHit().GetVertexIds(i) * 100000 + hit.GetTrueHit().GetPrimaryIds(i)] += hit.GetTrueHit().GetEnergyShare(i);
}
}
*/
} // End for (TG4HitSegmentContainer::iterator kt
} // End loop over each hit, for (TG4HitSegmentDetectors::iterator jt
}
Expand All @@ -282,15 +275,6 @@ TMS_Event::TMS_Event(TG4Event event, bool FillEvent) {
//std::cout<<"Making TMS event"<<std::endl;
bool OnlyPrimaryOrVisibleEnergy = true;

/*
if (LightWeight) {
OnlyMuon = true;
TMSOnly = true;
TMSLArOnly = true;
OnlyPrimary = true;
}
*/

// Save down the event number
EventNumber = EventCounter;
generator = std::default_random_engine(7890 + EventNumber);
Expand All @@ -305,12 +289,12 @@ TMS_Event::TMS_Event(TG4Event event, bool FillEvent) {

ProcessTG4Event(event, FillEvent);


// Now update truth info per particle
for (size_t i = 0; i < TMS_TrueParticles.size(); i++) {
double energy = 0;
// If it's not in the map, don't create it
auto it = TrueVisibleEnergyPerParticle.find(i);
int key = TMS_TrueParticles[i].GetVertexID() * 100000 + TMS_TrueParticles[i].GetTrackId();
auto it = TrueVisibleEnergyPerParticle.find(key);
if (it != TrueVisibleEnergyPerParticle.end()) {
energy = it->second;
}
Expand Down Expand Up @@ -388,14 +372,6 @@ void TMS_Event::MergeCoincidentHits() {
double y = (*it).GetNotZ();
//double e = hit.GetE();
double t = (*it).GetT();

// Strip out hits that are outside the actual volume
// This is probably some bug in the geometry that sometimes gives hits in the z=30k mm (i.e. 10m downstream of the end of the TMS)
// TODO figure out why these happen
if (z > TMS_Const::TMS_End[2] || z < TMS_Const::TMS_Start[2]) {
(*it).SetPedSup(true);
continue;
}

// Look ahead to find duplicates, but stop when z != z2
std::vector<std::vector<TMS_Hit>::iterator> duplicates;
Expand Down Expand Up @@ -886,6 +862,8 @@ void TMS_Event::AddEvent(TMS_Event &Other_Event) {
// Merge these lists
TrueVisibleEnergyPerVertex.merge(Other_Event.TrueVisibleEnergyPerVertex);
TrueVisibleEnergyPerParticle.merge(Other_Event.TrueVisibleEnergyPerParticle);
// Reset this to recalculate on next call
VertexIdOfMostEnergyInEvent = -9999;

nVertices += Other_Event.nVertices;

Expand Down Expand Up @@ -941,7 +919,7 @@ int TMS_Event::GetVertexIdOfMostVisibleEnergy() {
for (auto it : TrueVisibleEnergyPerVertex) {
double vertex = it.first;
double energy = it.second;
if (energy > max) { max = energy; max_vertex_id = vertex; }
if (energy >= max) { max = energy; max_vertex_id = vertex; }
total_energy += energy;
}
VertexIdOfMostEnergyInEvent = max_vertex_id;
Expand All @@ -950,7 +928,11 @@ int TMS_Event::GetVertexIdOfMostVisibleEnergy() {

if (TrueVisibleEnergyPerVertex.find(VertexIdOfMostEnergyInEvent) != TrueVisibleEnergyPerVertex.end())
TotalVisibleEnergyFromVertex = TrueVisibleEnergyPerVertex[VertexIdOfMostEnergyInEvent];
else std::cout<<"Warning in GetVertexIdOfMostVisibleEnergy: TrueVisibleEnergyPerVertex.Find(VertexIdOfMostEnergyInEvent) == TrueVisibleEnergyPerVertex.end()"<<std::endl;
else {
if (TrueVisibleEnergyPerVertex.size() > 0) {
std::cout<<"Warning in GetVertexIdOfMostVisibleEnergy: TrueVisibleEnergyPerVertex.Find(VertexIdOfMostEnergyInEvent) == TrueVisibleEnergyPerVertex.end()"<<std::endl;
}
}

return VertexIdOfMostEnergyInEvent;
}
Expand Down
Loading