Skip to content

[SC] Adding possibility to use float for calculations #6870

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

Merged
davidrohr merged 4 commits into from
Sep 1, 2021
Merged

[SC] Adding possibility to use float for calculations #6870

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f5fda7f
[SC] Adding possibility to use float for calculations
matthias-kleiner Aug 19, 2021
7a2e55a
fixing function name getChargeCyl -> getDensityCyl
matthias-kleiner Aug 19, 2021
43567bc
fixing typo
matthias-kleiner Aug 20, 2021
346ecd4
fixing error during compiling of macro
matthias-kleiner Aug 31, 2021
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: 2 additions & 2 deletions Detectors/TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -222,7 +222,7 @@ void debugInterpolation(utils::TreeStreamRedirector& pcstream,
float gxyzf[3] = {gx, gy, gz};
float pointCyl[3] = {gz, r, phi};
double efield[3] = {0.0, 0.0, 0.0};
double charge = spaceCharge->getChargeCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
double charge = spaceCharge->getDensityCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
double potential = spaceCharge->getPotentialCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
spaceCharge->getElectricFieldsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, efield[0], efield[1], efield[2]);
spaceCharge->getLocalDistortionsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, ldD[0], ldD[1], ldD[2]);
Expand Down Expand Up @@ -379,7 +379,7 @@ void debugGridpoints(utils::TreeStreamRedirector& pcstream, const o2::gpu::TPCFa
double efield[3] = {0.0, 0.0, 0.0};
double distLocal[3] = {0.0, 0.0, 0.0};
double dist[3] = {0.0, 0.0, 0.0};
double charge = spaceCharge->getChargeCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
double charge = spaceCharge->getDensityCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
double potential = spaceCharge->getPotentialCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
spaceCharge->getElectricFieldsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, efield[0], efield[1], efield[2]);
spaceCharge->getLocalDistortionsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, distLocal[0], distLocal[1], distLocal[2]);
Expand Down
26 changes: 13 additions & 13 deletions Detectors/TPC/spacecharge/include/TPCSpaceCharge/RegularGrid3D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,19 +137,19 @@ struct RegularGrid3D {

private:
inline static auto& mParamGrid = ParameterSpaceCharge::Instance();
static constexpr unsigned int FDIM = 3; ///< dimensions of the grid (only 3 supported)
static constexpr unsigned int FZ = 0; ///< index for x coordinate
static constexpr unsigned int FR = 1; ///< index for y coordinate
static constexpr unsigned int FPHI = 2; ///< index for z coordinate
const Vector<DataT, FDIM> mMin{}; ///< min vertices positions of the grid
const Vector<DataT, FDIM> mMax{}; ///< max vertices positions of the grid
const Vector<DataT, FDIM> mSpacing{}; ///< spacing of the grid
const Vector<DataT, FDIM> mInvSpacing{}; ///< inverse spacing of grid
const Vector<DataT, FDIM> sMaxIndex{{static_cast<double>(mParamGrid.NZVertices - 1.), static_cast<double>(mParamGrid.NRVertices - 1), static_cast<double>(mParamGrid.NPhiVertices - 1)}}; ///< max index which is on the grid in all dimensions
const Vector<int, FDIM> sNdim{{static_cast<int>(mParamGrid.NZVertices), static_cast<int>(mParamGrid.NRVertices), static_cast<int>(mParamGrid.NPhiVertices)}}; ///< number of vertices for each dimension
std::vector<DataT> mZVertices{}; ///< positions of vertices in x direction
std::vector<DataT> mRVertices{}; ///< positions of vertices in y direction
std::vector<DataT> mPhiVertices{}; ///< positions of vertices in z direction
static constexpr unsigned int FDIM = 3; ///< dimensions of the grid (only 3 supported)
static constexpr unsigned int FZ = 0; ///< index for x coordinate
static constexpr unsigned int FR = 1; ///< index for y coordinate
static constexpr unsigned int FPHI = 2; ///< index for z coordinate
const Vector<DataT, FDIM> mMin{}; ///< min vertices positions of the grid
const Vector<DataT, FDIM> mMax{}; ///< max vertices positions of the grid
const Vector<DataT, FDIM> mSpacing{}; ///< spacing of the grid
const Vector<DataT, FDIM> mInvSpacing{}; ///< inverse spacing of grid
const Vector<DataT, FDIM> sMaxIndex{{static_cast<DataT>(mParamGrid.NZVertices - 1.), static_cast<DataT>(mParamGrid.NRVertices - 1), static_cast<DataT>(mParamGrid.NPhiVertices - 1)}}; ///< max index which is on the grid in all dimensions
const Vector<int, FDIM> sNdim{{static_cast<int>(mParamGrid.NZVertices), static_cast<int>(mParamGrid.NRVertices), static_cast<int>(mParamGrid.NPhiVertices)}}; ///< number of vertices for each dimension
std::vector<DataT> mZVertices{}; ///< positions of vertices in x direction
std::vector<DataT> mRVertices{}; ///< positions of vertices in y direction
std::vector<DataT> mPhiVertices{}; ///< positions of vertices in z direction

void initLists();

Expand Down
29 changes: 24 additions & 5 deletions Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@

class TH3;
class TH3D;
class TH3F;

namespace o2
{
Expand All @@ -50,11 +51,15 @@ class SpaceCharge
using DataContainer = DataContainer3D<DataT>;
using GridProp = GridProperties<DataT>;
using TriCubic = TriCubicInterpolator<DataT>;
using TH3DataT = std::conditional_t<std::is_same<DataT, double>::value, TH3D, TH3F>; // datatype for TH3 (TH3F for DataT==float and TH3D for DataT==double)

public:
/// default constructor
/// grid granularity has to set before constructing an object using the static function setGrid(nZVertices, nRVertices, nPhiVertices)!
SpaceCharge() = default;
/// \param omegaTau \omega \tau value
/// \param t1 value for t1
/// \param t2 value for t2
SpaceCharge(const DataT omegaTau = -0.35f, const DataT t1 = 1, const DataT t2 = 1) { setOmegaTauT1T2(omegaTau, t1, t2); };

/// \param nZVertices number of vertices of the grid in z direction
/// \param nRVertices number of vertices of the grid in z direction
Expand Down Expand Up @@ -182,7 +187,7 @@ class SpaceCharge
/// \param z global z coordinate
/// \param r global r coordinate
/// \param phi global phi coordinate
DataT getChargeCyl(const DataT z, const DataT r, const DataT phi, const Side side) const;
DataT getDensityCyl(const DataT z, const DataT r, const DataT phi, const Side side) const;

/// get the potential for given coordinate
/// \param z global z coordinate
Expand Down Expand Up @@ -538,7 +543,7 @@ class SpaceCharge
/// \param omegaTau \omega \tau value
/// \param t1 value for t1 see: ???
/// \param t2 value for t2 see: ???
void setOmegaTauT1T2(const DataT omegaTau, const DataT t1, const DataT t2)
void setOmegaTauT1T2(const DataT omegaTau = -0.35f, const DataT t1 = 1, const DataT t2 = 1)
{
const DataT wt0 = t2 * omegaTau;
mC0 = 1 / (1 + wt0 * wt0);
Expand Down Expand Up @@ -567,7 +572,7 @@ class SpaceCharge
static void setNThreads(const int nThreads)
{
sNThreads = nThreads;
o2::tpc::TriCubicInterpolator<DataT>::setNThreads(nThreads);
TriCubic::setNThreads(nThreads);
}

/// set which kind of numerical integration is used for calcution of the integrals int Er/Ez dz, int Ephi/Ez dz, int Ez dz
Expand Down Expand Up @@ -597,6 +602,14 @@ class SpaceCharge
/// \side side of the TPC
void dumpToFile(TFile& outf, const Side side) const;

/// dump sc density, potential, electric fields, global distortions/corrections to tree
/// \param outFileName name of the output file
/// \side side of the TPC
/// \param nZPoints number of vertices of the output in z
/// \param nRPoints number of vertices of the output in r
/// \param nPhiPoints number of vertices of the output in phi
void dumpToTree(const char* outFileName, const Side side, const int nZPoints = 50, const int nRPoints = 150, const int nPhiPoints = 180) const;

/// write electric fields to root file
/// \param outf output file where the electrical fields will be written to
/// \side side of the TPC
Expand Down Expand Up @@ -716,7 +729,13 @@ class SpaceCharge
/// \param nBinsZNew number of z bins of rebinned histogram
/// \param nBinsRNew number of r bins of rebinned histogram
/// \param nBinsPhiNew number of phi bins of rebinned histogram
static TH3D rebinDensityHisto(const TH3& hOrig, const unsigned short nBinsZNew, const unsigned short nBinsRNew, const unsigned short nBinsPhiNew);
static TH3DataT rebinDensityHisto(const TH3& hOrig, const unsigned short nBinsZNew, const unsigned short nBinsRNew, const unsigned short nBinsPhiNew);

/// function to calculate the drift paths of the electron whose starting position is delivered. Electric fields must be set!
/// \param elePos global position of the start position of the electron
/// \param nSamplingPoints number of output points of the electron drift path
/// \return returns the input electron a vector of 3D-points describing the drift path of the electron
std::vector<GlobalPosition3D> calculateElectronDriftPath(const GlobalPosition3D& elePos, const int nSamplingPoints) const;

private:
inline static auto& mParamGrid = ParameterSpaceCharge::Instance(); ///< parameters of the grid on which the calculations are performed
Expand Down
3 changes: 2 additions & 1 deletion Detectors/TPC/spacecharge/include/TPCSpaceCharge/Vector.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,8 @@
#ifndef ALICEO2_TPC_VECTOR_H_
#define ALICEO2_TPC_VECTOR_H_

#include <Vc/Vc>
#include <Vc/vector>
#include <Vc/Memory>

namespace o2
{
Expand Down
48 changes: 32 additions & 16 deletions Detectors/TPC/spacecharge/macro/calculateDistortionsCorrections.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,11 +6,15 @@
#include "CommonUtils/TreeStreamRedirector.h"

template <typename DataT>
void calculateDistortionsAnalytical(const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1);
void calculateDistortionsAnalytical(const int sides = 0, const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1);

template <typename DataT>
void calculateDistortionsFromHist(const char* path, const char* histoName, const int sides, const int globalEFieldType = 1, const int globalDistType = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1);

int getSideStart(const int sides);
int getSideEnd(const int sides);

/// \param sides set which sides will be simulated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only
/// \param nZVertices number of vertices of the grid in z direction
/// \param nRVertices number of vertices of the grid in z direction
/// \param nPhiVertices number of vertices of the grid in z direction
Expand All @@ -21,10 +25,11 @@ void calculateDistortionsFromHist(const char* path, const char* histoName, const
/// \param nSteps number of which are used for calculation of distortions/corrections per z-bin
/// \param simpsonIterations number of iterations used in the simpson intergration
/// \param nThreads number of threads which are used (if the value is -1 all threads should be used)
void calcDistAna(const unsigned short nZVertices = 129, const unsigned short nRVertices = 129, const unsigned short nPhiVertices = 180, const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1)
template <typename DataT>
void calcDistAna(const int sides = 0, const unsigned short nZVertices = 129, const unsigned short nRVertices = 129, const unsigned short nPhiVertices = 180, const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1)
{
o2::tpc::SpaceCharge<double>::setGrid(nZVertices, nRVertices, nPhiVertices);
calculateDistortionsAnalytical<double>(globalEFieldTypeAna, globalDistTypeAna, eFieldTypeAna, usePoissonSolverAna, nSteps, simpsonIterations, nThreads);
o2::tpc::SpaceCharge<DataT>::setGrid(nZVertices, nRVertices, nPhiVertices);
calculateDistortionsAnalytical<DataT>(sides, globalEFieldTypeAna, globalDistTypeAna, eFieldTypeAna, usePoissonSolverAna, nSteps, simpsonIterations, nThreads);
}

/// \param path path to the root file containing the 3D density histogram
Expand Down Expand Up @@ -271,7 +276,7 @@ void writeToTree(o2::tpc::SpaceCharge<DataT>& spaceCharge3D, o2::utils::TreeStre
/// \param nSteps number of which are used for calculation of distortions/corrections per z-bin
/// \param simpsonIterations number of iterations used in the simpson intergration
template <typename DataT = double>
void calculateDistortionsAnalytical(const int globalEFieldTypeAna, const int globalDistTypeAna, const int eFieldTypeAna, const int usePoissonSolverAna, const int nSteps, const int simpsonIterations, const int nThreads)
void calculateDistortionsAnalytical(const int sides, const int globalEFieldTypeAna, const int globalDistTypeAna, const int eFieldTypeAna, const int usePoissonSolverAna, const int nSteps, const int simpsonIterations, const int nThreads)
{
const auto integrationStrategy = o2::tpc::SpaceCharge<DataT>::IntegrationStrategy::SimpsonIterative;
o2::tpc::SpaceCharge<DataT> spaceCharge3D;
Expand All @@ -286,7 +291,7 @@ void calculateDistortionsAnalytical(const int globalEFieldTypeAna, const int glo
// write to root file
o2::utils::TreeStreamRedirector pcstream(TString::Format("distortions_ana_nR%hu_nZ%hu_nPhi%hu_SimpsonsIter%i.root", spaceCharge3D.getNRVertices(), spaceCharge3D.getNZVertices(), spaceCharge3D.getNPhiVertices(), simpsonIterations).Data(), "RECREATE");

for (int iside = 0; iside < 2; ++iside) {
for (int iside = getSideStart(sides); iside < getSideEnd(sides); ++iside) {
std::cout << "side: " << iside << std::endl;
o2::tpc::Side side = (iside == 0) ? o2::tpc::Side::A : o2::tpc::Side::C;
o2::tpc::AnalyticalFields<DataT> anaFields(side);
Expand Down Expand Up @@ -327,16 +332,7 @@ void calculateDistortionsFromHist(const char* path, const char* histoName, const
fileOutSCDensity.Close();

o2::utils::TreeStreamRedirector pcstream(TString::Format("distortions_real_nR%hu_nZ%hu_nPhi%hu_SimpsonsIter%i.root", spaceCharge3D.getNRVertices(), spaceCharge3D.getNZVertices(), spaceCharge3D.getNPhiVertices(), simpsonIterations).Data(), "RECREATE");
int iSideStart = 0;
int iSideEnd = 2;
if (sides == 1) {
// a side only
iSideEnd = 1;
} else if (sides == 2) {
// c side only
iSideStart = 1;
}
for (int iside = iSideStart; iside < iSideEnd; ++iside) {
for (int iside = getSideStart(sides); iside < getSideEnd(sides); ++iside) {
std::cout << "side: " << iside << std::endl;
o2::tpc::Side side = (iside == 0) ? o2::tpc::Side::A : o2::tpc::Side::C;
const auto distType = globalDistType == 0 ? SC::GlobalDistType::Standard : SC ::GlobalDistType::Fast;
Expand Down Expand Up @@ -364,3 +360,23 @@ void calculateDistortionsFromHist(const char* path, const char* histoName, const
}
fOut.Close();
}

/// helper function to set the loop over the sides for the tpc
/// \param sides set for which sides the distortions/corrections will be calculated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only
int getSideStart(const int sides)
{
if (sides == 2) {
return 1;
}
return 0;
}

/// helper function to set the loop over the sides for the tpc
/// \param sides set for which sides the distortions/corrections will be calculated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only
int getSideEnd(const int sides)
{
if (sides == 1) {
return 1;
}
return 2;
}
Loading