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

TRestAxionMagneticField. nested std::vector to single std::vector approach #83

Closed
wants to merge 3 commits into from
Closed

TRestAxionMagneticField. nested std::vector to single std::vector approach #83

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
812b97c
TRestAxionMagneticField. Testing single std::vector performance
jgalan Feb 12, 2024
8e8727b
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Feb 12, 2024
10d21b3
Merge branch 'master' into field_upgrade
jgalan Feb 13, 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
4 changes: 3 additions & 1 deletion inc/TRestAxionMagneticField.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ struct MagneticFieldVolume {
TRestMesh mesh;

/// The field data connected to the grid defined by the mesh
std::vector<std::vector<std::vector<TVector3>>> field;
std::vector<TVector3> field;
};

/// A class to load magnetic field maps and evaluate the field on those maps including interpolation.
Expand Down Expand Up @@ -112,6 +112,8 @@ class TRestAxionMagneticField : public TRestMetadata {
}
}

size_t GetArrayElement(MagneticFieldVolume& mVol, Int_t nx, Int_t ny, Int_t nz);

public:
void LoadMagneticVolumes();

Expand Down
99 changes: 59 additions & 40 deletions src/TRestAxionMagneticField.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -696,6 +696,10 @@ TCanvas* TRestAxionMagneticField::DrawTracks(TVector3 vanishingPoint, Int_t divi
return fCanvas;
}

size_t TRestAxionMagneticField::GetArrayElement(MagneticFieldVolume& mVol, Int_t nx, Int_t ny, Int_t nz) {
return (size_t)(nx * mVol.mesh.GetNodesY() * mVol.mesh.GetNodesZ() + ny * mVol.mesh.GetNodesZ() + nz);
}

///////////////////////////////////////////////
/// \brief A method to help loading magnetic field data, as x,y,z,Bx,By,Bz into a magnetic volume
/// definition using its corresponding mesh.
Expand All @@ -704,12 +708,7 @@ TCanvas* TRestAxionMagneticField::DrawTracks(TVector3 vanishingPoint, Int_t divi
///
void TRestAxionMagneticField::LoadMagneticFieldData(MagneticFieldVolume& mVol,
std::vector<std::vector<Float_t>> data) {
mVol.field.resize(mVol.mesh.GetNodesX());
for (unsigned int n = 0; n < mVol.field.size(); n++) {
mVol.field[n].resize(mVol.mesh.GetNodesY());
for (unsigned int m = 0; m < mVol.field[n].size(); m++)
mVol.field[n][m].resize(mVol.mesh.GetNodesZ());
}
mVol.field.resize(mVol.mesh.GetNodesX() * mVol.mesh.GetNodesY() * mVol.mesh.GetNodesZ());

RESTDebug << "TRestAxionMagneticField::LoadMagneticFieldData. Printing first 5 data rows" << RESTendl;
for (unsigned int n = 0; n < data.size(); n++) {
Expand All @@ -720,6 +719,8 @@ void TRestAxionMagneticField::LoadMagneticFieldData(MagneticFieldVolume& mVol,
Int_t nY = mVol.mesh.GetNodeY((Int_t)(data[n][1] + mVol.mesh.GetNetSizeY() / 2.), true);
Int_t nZ = mVol.mesh.GetNodeZ((Int_t)(data[n][2] + mVol.mesh.GetNetSizeZ() / 2.), true);

size_t element = GetArrayElement(mVol, nX, nY, nZ);

if (n < 5) {
RESTDebug << "X: " << data[n][0] << " Y: " << data[n][1] << " Z: " << data[n][2] << RESTendl;
RESTDebug << "absX: " << data[n][0] + mVol.position.X()
Expand All @@ -730,15 +731,15 @@ void TRestAxionMagneticField::LoadMagneticFieldData(MagneticFieldVolume& mVol,
if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) GetChar();
}

if (mVol.field[nX][nY][nZ] != TVector3(0.0, 0.0, 0.0)) {
if (mVol.field[element] != TVector3(0.0, 0.0, 0.0)) {
RESTWarning << "X: " << data[n][0] << " Y: " << data[n][1] << " Z: " << data[n][2] << RESTendl;
RESTWarning << "nX: " << nX << " nY: " << nY << " nZ: " << nZ << RESTendl;
RESTWarning << "WARNING: field[nX][nY][nZ] element not equal to initial value (0, 0, 0) !!"
<< RESTendl;
RESTWarning << "It has value: "
<< "mVol.field[" << nX << "][" << nY << "][" << nZ << "] = ("
<< mVol.field[nX][nY][nZ].X() << " , " << mVol.field[nX][nY][nZ].Y() << " , "
<< mVol.field[nX][nY][nZ].Z() << ")" << RESTendl;
<< mVol.field[element].X() << " , " << mVol.field[element].Y() << " , "
<< mVol.field[element].Z() << ")" << RESTendl;
RESTWarning << "Values to write: "
<< "Bx: " << data[n][3] << " By: " << data[n][4] << " Bz: " << data[n][5] << RESTendl
<< RESTendl;
Expand All @@ -747,7 +748,7 @@ void TRestAxionMagneticField::LoadMagneticFieldData(MagneticFieldVolume& mVol,
if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) GetChar();
}

mVol.field[nX][nY][nZ] = TVector3(data[n][3], data[n][4], data[n][5]);
mVol.field[element] = TVector3(data[n][3], data[n][4], data[n][5]);
}
}

Expand Down Expand Up @@ -970,14 +971,30 @@ TVector3 TRestAxionMagneticField::GetMagneticField(TVector3 pos, Bool_t showWarn
else
nZ_1 = nZ;

TVector3 C000 = fMagneticFieldVolumes[id].field[nX][nY][nZ] + fConstantField[id];
TVector3 C100 = fMagneticFieldVolumes[id].field[nX_1][nY][nZ] + fConstantField[id];
TVector3 C010 = fMagneticFieldVolumes[id].field[nX][nY_1][nZ] + fConstantField[id];
TVector3 C110 = fMagneticFieldVolumes[id].field[nX_1][nY_1][nZ] + fConstantField[id];
TVector3 C001 = fMagneticFieldVolumes[id].field[nX][nY][nZ_1] + fConstantField[id];
TVector3 C101 = fMagneticFieldVolumes[id].field[nX_1][nY][nZ_1] + fConstantField[id];
TVector3 C011 = fMagneticFieldVolumes[id].field[nX][nY_1][nZ_1] + fConstantField[id];
TVector3 C111 = fMagneticFieldVolumes[id].field[nX_1][nY_1][nZ_1] + fConstantField[id];
size_t element;
element = GetArrayElement(fMagneticFieldVolumes[id], nX, nY, nZ);
TVector3 C000 = fMagneticFieldVolumes[id].field[element] + fConstantField[id];

element = GetArrayElement(fMagneticFieldVolumes[id], nX_1, nY, nZ);
TVector3 C100 = fMagneticFieldVolumes[id].field[element] + fConstantField[id];

element = GetArrayElement(fMagneticFieldVolumes[id], nX, nY_1, nZ);
TVector3 C010 = fMagneticFieldVolumes[id].field[element] + fConstantField[id];

element = GetArrayElement(fMagneticFieldVolumes[id], nX_1, nY_1, nZ);
TVector3 C110 = fMagneticFieldVolumes[id].field[element] + fConstantField[id];

element = GetArrayElement(fMagneticFieldVolumes[id], nX, nY, nZ_1);
TVector3 C001 = fMagneticFieldVolumes[id].field[element] + fConstantField[id];

element = GetArrayElement(fMagneticFieldVolumes[id], nX_1, nY, nZ_1);
TVector3 C101 = fMagneticFieldVolumes[id].field[element] + fConstantField[id];

element = GetArrayElement(fMagneticFieldVolumes[id], nX, nY_1, nZ_1);
TVector3 C011 = fMagneticFieldVolumes[id].field[element] + fConstantField[id];

element = GetArrayElement(fMagneticFieldVolumes[id], nX_1, nY_1, nZ_1);
TVector3 C111 = fMagneticFieldVolumes[id].field[element] + fConstantField[id];

Double_t x0 = fMagneticFieldVolumes[id].mesh.GetX(nX);
Double_t x1 = fMagneticFieldVolumes[id].mesh.GetX(nX_1);
Expand Down Expand Up @@ -1089,29 +1106,31 @@ void TRestAxionMagneticField::ReMap(const size_t& n, const TVector3& newMapSize)
return;
}

Int_t scaleX = (Int_t)(newMapSize.X() / fMeshSize[n].X());
Int_t scaleY = (Int_t)(newMapSize.Y() / fMeshSize[n].Y());
Int_t scaleZ = (Int_t)(newMapSize.Z() / fMeshSize[n].Z());

Int_t newNodesX = (fMagneticFieldVolumes[n].mesh.GetNodesX() - 1) / scaleX + 1;
Int_t newNodesY = (fMagneticFieldVolumes[n].mesh.GetNodesY() - 1) / scaleY + 1;
Int_t newNodesZ = (fMagneticFieldVolumes[n].mesh.GetNodesZ() - 1) / scaleZ + 1;

for (Int_t nx = 0; nx < newNodesX; nx++)
for (Int_t ny = 0; ny < newNodesY; ny++)
for (Int_t nz = 0; nz < newNodesZ; nz++)
fMagneticFieldVolumes[n].field[nx][ny][nz] =
fMagneticFieldVolumes[n].field[nx * scaleX][ny * scaleY][nz * scaleZ];

fMagneticFieldVolumes[n].mesh.SetNodes(newNodesX, newNodesY, newNodesZ);
fMagneticFieldVolumes[n].field.resize(fMagneticFieldVolumes[n].mesh.GetNodesX());
for (unsigned int i = 0; i < fMagneticFieldVolumes[n].field.size(); i++) {
fMagneticFieldVolumes[n].field[n].resize(fMagneticFieldVolumes[n].mesh.GetNodesY());
for (unsigned int j = 0; j < fMagneticFieldVolumes[n].field[i].size(); j++)
fMagneticFieldVolumes[n].field[i][j].resize(fMagneticFieldVolumes[n].mesh.GetNodesZ());
}
/*
Int_t scaleX = (Int_t)(newMapSize.X() / fMeshSize[n].X());
Int_t scaleY = (Int_t)(newMapSize.Y() / fMeshSize[n].Y());
Int_t scaleZ = (Int_t)(newMapSize.Z() / fMeshSize[n].Z());

Int_t newNodesX = (fMagneticFieldVolumes[n].mesh.GetNodesX() - 1) / scaleX + 1;
Int_t newNodesY = (fMagneticFieldVolumes[n].mesh.GetNodesY() - 1) / scaleY + 1;
Int_t newNodesZ = (fMagneticFieldVolumes[n].mesh.GetNodesZ() - 1) / scaleZ + 1;

for (Int_t nx = 0; nx < newNodesX; nx++)
for (Int_t ny = 0; ny < newNodesY; ny++)
for (Int_t nz = 0; nz < newNodesZ; nz++)
fMagneticFieldVolumes[n].field[nx][ny][nz] =
fMagneticFieldVolumes[n].field[nx * scaleX][ny * scaleY][nz * scaleZ];

fMagneticFieldVolumes[n].mesh.SetNodes(newNodesX, newNodesY, newNodesZ);
fMagneticFieldVolumes[n].field.resize(fMagneticFieldVolumes[n].mesh.GetNodesX());
for (unsigned int i = 0; i < fMagneticFieldVolumes[n].field.size(); i++) {
fMagneticFieldVolumes[n].field[n].resize(fMagneticFieldVolumes[n].mesh.GetNodesY());
for (unsigned int j = 0; j < fMagneticFieldVolumes[n].field[i].size(); j++)
fMagneticFieldVolumes[n].field[i][j].resize(fMagneticFieldVolumes[n].mesh.GetNodesZ());
}

fMeshSize[n] = TVector3(fMeshSize[n].X() * scaleX, fMeshSize[n].Y() * scaleY, fMeshSize[n].Z() * scaleZ);
fMeshSize[n] = TVector3(fMeshSize[n].X() * scaleX, fMeshSize[n].Y() * scaleY, fMeshSize[n].Z() * scaleZ);
*/
}

///////////////////////////////////////////////
Expand Down
Loading