Skip to content

Commit

Permalink
Merge pull request #30 from DUNE/kleykamp_larsoft_geometry_fix
Browse files Browse the repository at this point in the history 
Adding changes to scale geometry within dune-tms in the case of larsoft
  • Loading branch information
@LiamOS
LiamOS authored May 4, 2023
2 parents 29c1bc1 + 664bfe4 commit 656cf0f
Show file tree
Hide file tree
Showing 5 changed files with 149 additions and 30 deletions.
27 changes: 27 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
# Changelog

All notable changes to this project will be documented in this file.

The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
and this project adheres to a modified version of [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
Geometry releases will be tagged as `Descriptive_tag_v_X.Y.Z`.

## [Unreleased]
### Changed
- Fixes to make larsoft work. Specifically scaling.

## [Unreleased] - 2023-03-27
### Fixed
- Bugfix for infinite loop when a hit is outside volDetEnclosure

## [Unreleased] - 2023-03-17
### Added
- Timing simulation
- Time slicer
- Simple event display

## [TMSonlyFreeze] - 2022-03-01
- Tag for frozen TMS reconstruction. Used for the flux studies on TMS only

## [TrackMatching] - 2021-06-28
- Used for track matching studies with third production
8 changes: 5 additions & 3 deletions src/TMS_Bar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ bool TMS_Bar::FindModules(double xval, double yval, double zval) {
TGeoManager *geom = TMS_Geom::GetInstance().GetGeometry();

// Find which node this position is equivalent too
std::string NodeName = std::string(geom->FindNode(xval,yval,zval)->GetName());
std::string NodeName = std::string(TMS_Geom::GetInstance().FindNode(xval,yval,zval)->GetName());

// The position of the hit bar
double Translation[] = {0., 0., 0.};
Expand All @@ -57,6 +57,7 @@ bool TMS_Bar::FindModules(double xval, double yval, double zval) {
xw = 2*box->GetDX();
yw = 2*box->GetDY();
zw = 2*box->GetDZ();
TMS_Geom::GetInstance().Scale(xw, yw, zw);

// Do a sanity check (CHEATING!)
// Know the bars are 1cm in z and 3.542cm in x
Expand All @@ -68,6 +69,7 @@ bool TMS_Bar::FindModules(double xval, double yval, double zval) {
}

for (int i = 0; i < 3; ++i) Translation[i] = geom->GetCurrentMatrix()->GetTranslation()[i];
TMS_Geom::GetInstance().Scale(Translation[0], Translation[1], Translation[2]);
}

else if (NodeName.find(TMS_Const::TMS_ModuleName) != std::string::npos) {
Expand Down Expand Up @@ -118,7 +120,7 @@ bool TMS_Bar::FindModules(double xval, double yval, double zval) {
}

// Reset the geom navigator node level in case it's used again
geom->FindNode(xval,yval,zval);
TMS_Geom::GetInstance().FindNode(xval,yval,zval);

// Update the bar number
BarNumber = (GetX()-TMS_Const::TMS_Start_Exact[0])/GetXw();
Expand All @@ -135,7 +137,7 @@ int TMS_Bar::FindBar(double x, double y, double z) {
TGeoManager *geom = TMS_Geom::GetInstance().GetGeometry();

// Find which node this position is equivalent too
geom->FindNode(x,y,z);
TMS_Geom::GetInstance().FindNode(x,y,z);
TGeoNavigator *nav = geom->GetCurrentNavigator();
std::string NodeName = std::string(nav->GetCurrentNode()->GetName());
// cd up in the geometry to find the right name
Expand Down
9 changes: 8 additions & 1 deletion src/TMS_Event.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) {
TG4TrajectoryPoint pt = *kt;

// Check the point against the geometry
TGeoNode *vol = TMS_Geom::GetInstance().GetGeometry()->FindNode(pt.GetPosition().X(), pt.GetPosition().Y(), pt.GetPosition().Z());
TGeoNode *vol = TMS_Geom::GetInstance().FindNode(pt.GetPosition().X(), pt.GetPosition().Y(), pt.GetPosition().Z());

// Very rarely but it does happen, the volume is null
if (!vol) continue;
Expand Down Expand Up @@ -212,6 +212,7 @@ TMS_Event::TMS_Event(TG4Event &event, bool FillEvent) {
} // End for (TG4HitSegmentContainer::iterator kt
} // End loop over each hit, for (TG4HitSegmentDetectors::iterator jt


// Merge hits that happened in the same scintillator strip and within the same readout time window
// This is a simulation cleanup step, not reconstruction
MergeCoincidentHits();
Expand Down Expand Up @@ -292,6 +293,7 @@ void TMS_Event::MergeCoincidentHits() {
for (auto hit : TMS_Hits) {
if (!hit.GetPedSup()) remaining_hits.push_back(hit);
else deleted_hits.push_back(hit);
if (!hit.GetPedSup() && hit.GetE() > 10000) std::cout << "Warning: Found hit higher than 10 GeV. Seems unlikely. Hit E = " << (hit.GetE() / 1000.0) << " GeV." << std::endl;
}
TMS_Hits.clear();
for (auto hit : remaining_hits) TMS_Hits.push_back(hit);
Expand Down Expand Up @@ -357,6 +359,11 @@ void TMS_Event::SimulateTimingModel() {
double time_correction_long_way = long_way_distance / SPEED_OF_LIGHT_IN_FIBER;
// Time slew (up to 30ns for 1pe hits, 9ns for 5pe, ~2ns 22pe. Typically 22pe mips assuming 45 pe mips with half going the long way)
double pe = hit.GetPE();
// We don't have to do 1000s of throws. The time will be very close to zero.
// Assuming 1k PE, the mean time is ~0.02ns vs ~0.06ns for 300 PE.
if (pe > 300) {
pe = 300;
}
double minimum_time_offset = 1e100;
while (pe > 0) {
// Light can either go the directly to the readout or go the long way first.
Expand Down
125 changes: 101 additions & 24 deletions src/TMS_Geom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#include "TGeoManager.h"

#include "TVector3.h"
#include "TGeoBBox.h"

#include "CLHEP/Units/SystemOfUnits.h"

Expand Down Expand Up @@ -46,16 +47,76 @@ class TMS_Geom {
// Set the geometry
void SetGeometry(TGeoManager *geometry) {
geom = geometry;
std::cout << "Global geometry set to " << geometry->GetName() << std::endl;
geom->LockGeometry();

// There's an overall scale factor depending on if the gmdl was loaded with cm or mm.
// So here we're trying to automatically figure out the scale factor
// If the geometry changes too much, then this would not work and so we'd want to give up.
TGeoBBox *box = dynamic_cast<TGeoBBox*>(geom->GetTopVolume()->GetShape());
double dx = 2*box->GetDX();
if (dx == 600000) ScaleFactor = 1;
else if (dx == 60000) ScaleFactor = 10;
else {
std::cerr << "Fatal: Unable to guess geometry's scale factor based on Shape for geometry " << geometry->GetName() << std::endl;
throw;
}
std::cout << "Global geometry set to " << geometry->GetName() << std::endl;
std::cout << "Geometry scale factor: " << ScaleFactor;
std::cout << ". Factor is 1 if 1 unit = 1mm (aka edep sim), 10 if 1 unit = 1cm (aka larsoft)." << std::endl;
}

void SetFileName(std::string filename) {
FileName = filename;
}

inline void Scale(double &x, double &y, double &z) {
x *= ScaleFactor;
y *= ScaleFactor;
z *= ScaleFactor;
}

inline double Scale(double val) {
return val * ScaleFactor;
}

TVector3 Scale(TVector3 point) {
double x = point.X();
double y = point.Y();
double z = point.Z();
Scale(x, y, z);
TVector3 out(x, y, z);
return out;
}

inline void Unscale(double &x, double &y, double &z) {
x /= ScaleFactor;
y /= ScaleFactor;
z /= ScaleFactor;
}

inline double Unscale(double val) {
return val / ScaleFactor;
}

TVector3 Unscale(TVector3 point) {
double x = point.X();
double y = point.Y();
double z = point.Z();
Unscale(x, y, z);
TVector3 out(x, y, z);
return out;
}

TGeoNode* FindNode(double x, double y, double z) {
Unscale(x, y, z);
return geom->FindNode(x, y, z);
}

// Largely modlled on TGeoChecker::ShootRay in ROOT (except there's a stopping point, not an infinitely long ray)
std::vector<std::pair<TGeoMaterial*, double> > GetMaterials(const TVector3 &point1, const TVector3 &point2) {
std::vector<std::pair<TGeoMaterial*, double> > GetMaterials(const TVector3 &point1_temp, const TVector3 &point2_temp) {
// Make vectors have geometry scale
TVector3 point1 = Unscale(point1_temp);
TVector3 point2 = Unscale(point2_temp);

// First cd the navigator to the starting point
geom->FindNode(point1.X(), point1.Y(), point1.Z());
Expand All @@ -82,7 +143,7 @@ class TMS_Geom {
// Count up the total length for debugging
double total = 0;
// Walk through until we're in the same volume as our final point
while (step < __GEOM_LARGE_STEP__ && target_dist-dist > 0) {
while (step < Unscale(__GEOM_LARGE_STEP__) && target_dist-dist > 0) {
// Get the material of the current point
TGeoMaterial *mat = geom->GetCurrentNode()->GetMedium()->GetMaterial();
// Get the position
Expand All @@ -98,8 +159,8 @@ class TMS_Geom {
geom->FindNode();
// How big was the step in this material
step = geom->GetStep();
if (step < __GEOM_TINY_STEP__) {
geom->SetStep(__GEOM_SMALL_STEP__);
if (step < Unscale(__GEOM_TINY_STEP__)) {
geom->SetStep(Unscale(__GEOM_SMALL_STEP__));
// Step into the next volume
geom->Step();
// Go down to the deepest node
Expand All @@ -114,7 +175,7 @@ class TMS_Geom {
}

// Push back the information
std::pair<TGeoMaterial*, double> temp(mat, step);
std::pair<TGeoMaterial*, double> temp(mat, Scale(step));
Materials.push_back(temp);
total += step;
}
Expand Down Expand Up @@ -154,7 +215,11 @@ class TMS_Geom {
return Materials;
}

std::vector<std::pair<TGeoMaterial*, TVector3> > GetMaterialsPos(const TVector3 &point1, const TVector3 &point2) {
std::vector<std::pair<TGeoMaterial*, TVector3> > GetMaterialsPos(const TVector3 &point1_temp, const TVector3 &point2_temp) {
// Make vectors have geometry scale
TVector3 point1 = Unscale(point1_temp);
TVector3 point2 = Unscale(point2_temp);

// First cd the navigator to the starting point
geom->FindNode(point1.X(), point1.Y(), point1.Z());

Expand All @@ -177,8 +242,8 @@ class TMS_Geom {

double step = geom->GetStep();
// Walk through until we're in the same volume as our final point
//while (!geom->IsSameLocation(point2.X(), point2.Y(), point2.Z()) && step < __GEOM_LARGE_STEP__) {
while (step < __GEOM_LARGE_STEP__ && target_dist-dist > 0) {
//while (!geom->IsSameLocation(point2.X(), point2.Y(), point2.Z()) && step < Unscale(__GEOM_LARGE_STEP__)) {
while (step < Unscale(__GEOM_LARGE_STEP__) && target_dist-dist > 0) {
// Get the material of the current point
TGeoMaterial *mat = geom->GetCurrentNode()->GetMedium()->GetMaterial();
// Get the position
Expand All @@ -195,8 +260,8 @@ class TMS_Geom {
// Check IsEntering?
// If IsEntering, set step to very small until not entering anymore
step = geom->GetStep();
if (step < __GEOM_TINY_STEP__) {
geom->SetStep(__GEOM_SMALL_STEP__);
if (step < Unscale(__GEOM_TINY_STEP__)) {
geom->SetStep(Unscale(__GEOM_SMALL_STEP__));
// Step into the next volume
geom->Step();
// Go down to the deepest node
Expand All @@ -208,14 +273,18 @@ class TMS_Geom {
if (dist+step > target_dist) break;

// Push back the information
std::pair<TGeoMaterial*, TVector3> temp(mat, pt_vec);
std::pair<TGeoMaterial*, TVector3> temp(mat, Scale(pt_vec));
Materials.push_back(temp);
}
return Materials;
}

// Get the track length between two points by walking through the materials at those points
double GetTrackLength(const TVector3 &point1, const TVector3 &point2) {
double GetTrackLength(const TVector3 &point1_temp, const TVector3 &point2_temp) {
// Make vectors have geometry scale
TVector3 point1 = Unscale(point1_temp);
TVector3 point2 = Unscale(point2_temp);

// First get the collection of materials between point1 and point2
std::vector<std::pair<TGeoMaterial*, double> > Materials = GetMaterials(point1, point2);

Expand Down Expand Up @@ -246,11 +315,14 @@ class TMS_Geom {

counter++;
}

TotalPathLength = Scale(TotalPathLength);
return TotalPathLength;
};

std::vector<std::pair<std::string, TVector3> > GetNodes(const TVector3 &point1, const TVector3 &point2) {
std::vector<std::pair<std::string, TVector3> > GetNodes(const TVector3 &point1_temp, const TVector3 &point2_temp) {
// Make vectors have geometry scale
TVector3 point1 = Unscale(point1_temp);
TVector3 point2 = Unscale(point2_temp);

// First cd the navigator to the starting point
geom->FindNode(point1.X(), point1.Y(), point1.Z());
Expand All @@ -273,7 +345,7 @@ class TMS_Geom {
double step = geom->GetStep();

// Walk through until we're in the same volume as our final point
while (step < __GEOM_LARGE_STEP__ && target_dist-dist > 0) {
while (step < Unscale(__GEOM_LARGE_STEP__) && target_dist-dist > 0) {
// Get the material of the current point
std::string nodename = std::string(geom->GetCurrentNode()->GetName());
const double *pt = geom->GetCurrentPoint();
Expand All @@ -294,8 +366,8 @@ class TMS_Geom {
std::cout << geom->GetCurrentNode()->GetName() << std::endl;
// How big was the step in this material
step = geom->GetStep();
if (step < __GEOM_TINY_STEP__) {
geom->SetStep(__GEOM_SMALL_STEP__);
if (step < Unscale(__GEOM_TINY_STEP__)) {
geom->SetStep(Unscale(__GEOM_SMALL_STEP__));
// Step into the next volume
geom->Step();
// Go down to the deepest node
Expand All @@ -312,13 +384,16 @@ class TMS_Geom {
}

// Push back the information
std::pair<std::string, TVector3> temp(nodename, pt_vec);
std::pair<std::string, TVector3> temp(nodename, Scale(pt_vec));
Nodes.push_back(temp);
}
return Nodes;
}

std::vector<std::pair<int*, TVector3> > GetUniquePlaneBarIdent(const TVector3 &point1, const TVector3 &point2) {
std::vector<std::pair<int*, TVector3> > GetUniquePlaneBarIdent(const TVector3 &point1_temp, const TVector3 &point2_temp) {
// Make vectors have geometry scale
TVector3 point1 = Unscale(point1_temp);
TVector3 point2 = Unscale(point2_temp);

// First cd the navigator to the starting point
geom->FindNode(point1.X(), point1.Y(), point1.Z());
Expand All @@ -341,7 +416,7 @@ class TMS_Geom {
double step = geom->GetStep();

// Walk through until we're in the same volume as our final point
while (step < __GEOM_LARGE_STEP__ && target_dist-dist > 0) {
while (step < Unscale(__GEOM_LARGE_STEP__) && target_dist-dist > 0) {
// Plane, bar, global
int *Plane = new int[3];
for (int i = 0; i < 3; ++i) Plane[i] = -999;
Expand Down Expand Up @@ -375,8 +450,8 @@ class TMS_Geom {
geom->FindNextBoundaryAndStep();
// How big was the step in this material
step = geom->GetStep();
if (step < __GEOM_TINY_STEP__) {
geom->SetStep(__GEOM_SMALL_STEP__);
if (step < Unscale(__GEOM_TINY_STEP__)) {
geom->SetStep(Unscale(__GEOM_SMALL_STEP__));
// Step into the next volume
geom->Step();
// Go down to the deepest node
Expand All @@ -389,7 +464,7 @@ class TMS_Geom {
if (Plane[0] == -999 || Plane[1] == -999 || Plane[2] == -999) continue;

// Push back the information
std::pair<int*, TVector3> temp(Plane, pt_vec);
std::pair<int*, TVector3> temp(Plane, Scale(pt_vec));
Nodes.push_back(temp);
}

Expand All @@ -402,6 +477,7 @@ class TMS_Geom {
TMS_Geom() {
FileName = TMS_Manager::GetInstance().GetFileName();
geom = NULL;
ScaleFactor = 1;
};

~TMS_Geom() {};
Expand All @@ -412,6 +488,7 @@ class TMS_Geom {
// The actual geometry
TGeoManager *geom;
std::string FileName;
double ScaleFactor;
};

#endif
10 changes: 8 additions & 2 deletions src/TMS_Kalman.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,8 +194,14 @@ void TMS_Kalman::Predict(TMS_KalmanNode &Node) {
for (auto material : Materials) {

// Read these directly from a TGeoManager
double density = material.first->GetDensity()/(CLHEP::g/CLHEP::cm3); // now in g/cm3 (edep-sim geometry is in CLHEP units)
double thickness = material.second/10.; // in cm (was in mm in geometry)
// If the geometry is in mm (CLHEP units), then want to scale density to g/cm3 and thickness to cm
// Otherwise, assume it's like that and then fix it with geometry scaling functions
double density = material.first->GetDensity()/(CLHEP::g/CLHEP::cm3);
double thickness = material.second/10.;
// Potentially need to scale from g/cm3 to g/mm3, so find the scale factor and scale by 1/that^3.
double scale_factor = TMS_Geom::GetInstance().Scale(1.0);
density /= std::pow(scale_factor, 3);
thickness = TMS_Geom::GetInstance().Scale(thickness);
material.first->Print();

if (Talk) {
Expand Down

0 comments on commit 656cf0f

Please sign in to comment.