Adding classes to calculate the space charge distortions and corrections #4618

matthias-kleiner · 2020-10-16T17:09:03Z

The distortions and corrections can be calculated by using the potential
from a 3D histogram or an analytical formula as an input. These classes are
much faster and consistent than the current methods located in GPU/TPCSpaceChargeBase and will replace them.

wiechula

Perhaps instead of ...Structs.h you could use ...Helpers.h. In some of the file you declare structs with private members. I think I would rather make them classes.

Otherwise only minor comments

wiechula · 2020-10-20T07:34:56Z

Detectors/TPC/spacecharge/include/TPCSpacecharge/PoissonSolver.h

+  /// Helper function to check if the integer is equal to a power of two
+  /// \param i the number
+  /// \return 1 if it is a power of two, else 0
+  bool isPowerOfTwo(int i) const;


This could be simply

return ((i>0) && !(i & (i - 1)));

wiechula · 2020-10-20T07:36:44Z

Detectors/TPC/spacecharge/include/TPCSpacecharge/PoissonSolverStructs.h

+  /// operator to set the values
+  DataT& operator()(const unsigned int iR, const unsigned int iZ, const unsigned int iPhi)
+  {
+    return mStorage[iR + mNr * (iZ + mNz * iPhi)];


why not use getIndex

wiechula · 2020-10-20T07:36:49Z

Detectors/TPC/spacecharge/include/TPCSpacecharge/PoissonSolverStructs.h

+  /// operator to read the values
+  const DataT& operator()(const unsigned int iR, const unsigned int iZ, const unsigned int iPhi) const
+  {
+    return mStorage[iR + mNr * (iZ + mNz * iPhi)];


why not use getIndex

wiechula · 2020-10-20T07:38:26Z

Detectors/TPC/spacecharge/include/TPCSpacecharge/PoissonSolverStructs.h

+/// this is a simple vector class which is used in the poisson solver class
+
+/// \tparam DataT the data type of the mStorage which is used during the calculations
+template <typename DataT = double>


Why is this not in an own header as e.g. DataContainer3D?

wiechula · 2020-10-20T07:50:26Z

Detectors/TPC/spacecharge/include/TPCSpacecharge/SpaceChargeStructs.h

+#include <cmath>
+#include "TPCSpacecharge/TriCubic.h"
+#include "DataFormatsTPC/Defs.h"
+#include "TPCSpacecharge/SpaceChargeStructs.h"


wiechula · 2020-10-20T08:00:49Z

Detectors/TPC/spacecharge/include/TPCSpacecharge/SpaceChargeStructs.h

+
+  static constexpr unsigned int getID() { return ID; }
+
+  // private:


why commented?

wiechula · 2020-10-20T08:10:06Z

Detectors/TPC/spacecharge/include/TPCSpacecharge/Vector.h

+
+ public:
+  /// default constructor
+  Vector() {}


wiechula · 2020-10-20T08:15:54Z

Detectors/TPC/spacecharge/src/SpaceCharge.cxx

+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+  typedef Kernel::Point_3 Point_3;
+  typedef boost::tuple<Point_3, std::array<unsigned int, 3>> Point_Index;
+  typedef CGAL::Search_traits_3<Kernel> Traits_base;
+  typedef CGAL::Search_traits_adapter<Point_Index, CGAL::Nth_of_tuple_property_map<0, Point_Index>, Traits_base> Traits;
+  typedef CGAL::Orthogonal_k_neighbor_search<Traits> K_neighbor_search;
+  typedef K_neighbor_search::Tree Tree;


change to using?

ehellbar · 2020-10-21T14:11:50Z

Detectors/TPC/spacecharge/include/TPCSpacecharge/PoissonSolverStructs.h

+
+template <typename DataT = double>
+struct TPCParameters {
+  static constexpr DataT TPCZ0{249.525};                        ///< nominal gating grid position


should be nominal G1T position in the comment

ehellbar · 2020-10-21T14:29:52Z

Detectors/TPC/spacecharge/src/PoissonSolver.cxx

+void PoissonSolver<DataT, Nz, Nr, Nphi>::poissonSolver2D(DataContainer& matricesV, const DataContainer& matricesCharge)
+{
+  poissonMultiGrid2D(matricesV, matricesCharge);
+}


Just a note:
iPhi (third argument of the function) is set to 0 by default. In the function itself, only the iPhi bin of the output matrix matricesV will be filled (line 181). So for the 2D solution, only 1 phi bin is calculated and filled. Once (if) we start to use the 2D solutions for the ion drift, we have to revisit this.

ehellbar · 2020-10-21T16:23:59Z

Detectors/TPC/spacecharge/test/testO2TPCPoissonSolver.cxx

+  const o2::tpc::AnalyticalFields<DataT> analyticalFields;
+  // set the boudnary and charge for numerical poisson solver
+  setChargeDensityFromFormula<DataT, Nz, Nr, Nphi>(analyticalFields, grid3D, charge);
+  setPotentialBoundaryFromFormula<DataT, Nz, Nr, Nphi>(analyticalFields, grid3D, potentialNumerical);


Is there a reason why you set the boundary potentials of the numerical calculation with the values from the analytical formula? The numerical solution should work fine without this as it should calculate the potential at these boundaries from the charge?

The function to set only the boundary potentials is supposed to provide to possibility to account for irregularities of the nominal potentials, e.g. charge-up effects at dielectric sector boundaries.

Is there a reason why you set the boundary potentials of the numerical calculation with the values from the analytical formula?

It was like this in the original code. So I always used theses as the boundary conditions.

The numerical solution should work fine without this as it should calculate the potential at these boundaries from the charge?

I am not sure if boundary conditions are really needed or not for this method. I made a quick comparison by setting the boundary potential and then calculating the potential with the poisonsolver and without setting the boundary potential. The difference is quit big between the resulting potential to the reference potential. See the two plots "analytical potential - numerical potential". The first plot is with the boundary conditions and the setting plot without these boundary conditions for the potential.

Ok. Now I am thinking about what this means for the case when we use density histograms as input. Maybe at some point we can try to put there the nominal TPC boundary conditions (potential CE, G1T and field cage) already for the Poisson solver instead only adding it (Ez0=400 V/cm) for the calculation of the distortions.

ehellbar · 2020-10-21T16:34:56Z

Hello @matthias-kleiner ,

can you add a method to the main class which calculates all the quantities (potential, E-field, distortion, correction), similar to the function calculateDistortionsCorrections of the example macro? It would make the application of the class in future studies much simpler, or what do you think?

Otherwise, it looks very good to me.

Cheers,
Ernst

matthias-kleiner · 2020-10-23T15:35:26Z

Hello @matthias-kleiner ,

can you add a method to the main class which calculates all the quantities (potential, E-field, distortion, correction), similar to the function calculateDistortionsCorrections of the example macro? It would make the application of the class in future studies much simpler, or what do you think?

Otherwise, it looks very good to me.

Cheers,
Ernst

Hi Ernst,
I just implemented such a method and use this method now in the macro.

I also implemented all the other comments from you and Jens.

matthias-kleiner · 2020-10-23T15:50:15Z

@ehellbar

I still have to invert the distortions in z direction for the A-(?)Side. Should I do it directly during the filling of the values or when quering the values (like it was done in the old code). I assume it should be done during the filling?

And there is in the TPCSimulation folder a class with the same name "SpaceCharge.h", which causes some problems when running a simulation.

ehellbar · 2020-10-25T16:46:27Z

Hi @matthias-kleiner ,

I still have to invert the distortions in z direction for the A-(?)Side. Should I do it directly during the filling of the values or when quering the values (like it was done in the old code). I assume it should be done during the filling?

I think it concerns the values on the C side? On the A side, the direction along z in the array/matrix entries is 0 -> nZ and for the z coordinate 0 -> 250 cm. On the C side it is 0 -> nZ and 0 -> -250 cm. Differentiation / integration along z is done using the array / matrix index, right? So on the C side, you will basically invert the z axis (direction -z instead of z) which results in a change of the sign. I think this happens first when calculating the Ez field from the potential? So you can change the sign of Ez on the C side and it should be correct for all following steps? But we should check the signs of the E fields, local and global distortions after the change.

And there is in the TPCSimulation folder a class with the same name "SpaceCharge.h", which causes some problems when running a simulation.

Is it related to the new classes? They are not yet used in the simulation and are in a different scope, right (e.g. TPCSpaceCharge instead of TPCSimulation)?

Cheers,
Ernst

matthias-kleiner · 2020-10-26T12:54:11Z

Hi @ehellbar

Is it related to the new classes? They are not yet used in the simulation and are in a different scope, right (e.g. TPCSpaceCharge instead of TPCSimulation)?

when running a simulation it completes, but I get some error from I think ROOT:
error: redefinition of 'SpaceCharge' as different kind of symbol
note: previous definition is here
namespace o2{namespace tpc{class __attribute__((annotate("$clingAutoload$TPCSimulation/SpaceCharge.h"))) SpaceCharge;}}
I think its because they are in the same namespace o2::tpc

matthias-kleiner · 2020-10-26T14:16:43Z

Hi @matthias-kleiner ,

I still have to invert the distortions in z direction for the A-(?)Side. Should I do it directly during the filling of the values or when quering the values (like it was done in the old code). I assume it should be done during the filling?

I think it concerns the values on the C side?

Hi @ehellbar,
you are right. Its the value on the C-Side.

So you can change the sign of Ez on the C side and it should be correct for all following steps? But we should check the signs of the E fields, local and global distortions after the change.

I just checked it, but the resulting corrections/distortions are slightly different after changing the sign of the ez field. I think changing the sign of Ez also changes these integrals "int er/ez dz" etc. when calculating the distortions. So i think one should integrate then from 0 -> -nZBins, but then I have to change several other things too.

ehellbar · 2020-10-26T15:14:33Z

Hello @matthias-kleiner ,
yes, you are right. You have to check the sign everywhere where the Ez is used. In the integrals Er/(ezField+Ez), Ephi/(ezField+Ez) one has to change also the sign of ezField accordingly, as it has negative sign on the A side and positive sign on the C side (?).

ehellbar · 2020-10-26T15:42:51Z

when running a simulation it completes, but I get some error from I think ROOT:
error: redefinition of 'SpaceCharge' as different kind of symbol
note: previous definition is here
namespace o2{namespace tpc{class attribute((annotate("$clingAutoload$TPCSimulation/SpaceCharge.h"))) SpaceCharge;}}
I think its because they are in the same namespace o2::tpc

The class is mainly used in the TPC digitizer and some macros.

(base) hellbaer@hellbaer-Thinkpad:~/software/alice/o2-dev/O2$ grep -nr "TPCSimulation/SpaceCharge" *
Detectors/TPC/simulation/src/SpaceCharge.cxx:28:#include "TPCSimulation/SpaceCharge.h"
Detectors/TPC/simulation/CMakeLists.txt:42:                                  include/TPCSimulation/SpaceCharge.h)
Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h:21:#include "TPCSimulation/SpaceCharge.h"
Detectors/TPC/simulation/macro/createResidualDistortionObject.C:25:#include "TPCSimulation/SpaceCharge.h"
GPU/TPCFastTransformation/macro/generateTPCCorrectionNTuple.C:27:#include "TPCSimulation/SpaceCharge.h"

Maybe the best solution at the moment is to replace the usage of the other SpaceCharge class by your code, which requires a couple of new methods used in the other SpaceCharge class: distortElectron and correctElectron and an enum class SCDistortionType. Then, the old class can be removed from compilation. It still contains methods to propagate the space-charge density or ions which we should continue to develop and move to the spacecharge folder, either as a new class or as addition to SpaceCharge. But once the dependencies on the old SpaceCharge class are replaced by yours, we can basically delete the old one (getting the ion movement methods later from the history) and the GPU/TPCSpaceChargeBase folder as well.

matthias-kleiner · 2020-10-30T10:52:32Z

when running a simulation it completes, but I get some error from I think ROOT:
error: redefinition of 'SpaceCharge' as different kind of symbol
note: previous definition is here
namespace o2{namespace tpc{class attribute((annotate("$clingAutoload$TPCSimulation/SpaceCharge.h"))) SpaceCharge;}}
I think its because they are in the same namespace o2::tpc

The class is mainly used in the TPC digitizer and some macros.
(base) hellbaer@hellbaer-Thinkpad:~/software/alice/o2-dev/O2$ grep -nr "TPCSimulation/SpaceCharge" *
Detectors/TPC/simulation/src/SpaceCharge.cxx:28:#include "TPCSimulation/SpaceCharge.h"
Detectors/TPC/simulation/CMakeLists.txt:42:                                  include/TPCSimulation/SpaceCharge.h)
Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h:21:#include "TPCSimulation/SpaceCharge.h"
Detectors/TPC/simulation/macro/createResidualDistortionObject.C:25:#include "TPCSimulation/SpaceCharge.h"
GPU/TPCFastTransformation/macro/generateTPCCorrectionNTuple.C:27:#include "TPCSimulation/SpaceCharge.h"
Maybe the best solution at the moment is to replace the usage of the other SpaceCharge class by your code, which requires a couple of new methods used in the other SpaceCharge class: distortElectron and correctElectron and an enum class SCDistortionType. Then, the old class can be removed from compilation. It still contains methods to propagate the space-charge density or ions which we should continue to develop and move to the spacecharge folder, either as a new class or as addition to SpaceCharge. But once the dependencies on the old SpaceCharge class are replaced by yours, we can basically delete the old one (getting the ion movement methods later from the history) and the GPU/TPCSpaceChargeBase folder as well.

Hi @ehellbar,
since the new SpaceCharge class is templated, I am not sure about the best way to pass the template parameters when one uses this class in the Digitizer. Of course I can set a default grid size, but then one looses the flexibility to change the gridsize during run time. I could just get rid of all the templates and use vectors everywhere, but that would probably require some modifications of the code. What do you think about it?

ehellbar · 2020-10-30T15:19:51Z

Hi @matthias-kleiner ,
for the current use in the digitizer we can simply fix the grid size, e.g. to 129 x 129 x 180 (r x z x phi). The calculations with this grid size are quite fast (~10 sec?) and should be sufficiently precise. For studies e.g. of different grid sizes, we should use the class standalone. Fixing the grid size requires to also remove the option "gridSize" in the digitizer workflow O2/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx . Also, we have to check and modify the options initialSpaceChargeDensity (providing file name and histogram name) to calculate the distortion from a histogram and readSpaceCharge (provide a precalculated SpaceCharge object).

matthias-kleiner · 2020-10-30T17:20:00Z

Hi @matthias-kleiner ,
for the current use in the digitizer we can simply fix the grid size, e.g. to 129 x 129 x 180 (r x z x phi). The calculations with this grid size are quite fast (~10 sec?) and should be sufficiently precise. For studies e.g. of different grid sizes, we should use the class standalone. Fixing the grid size requires to also remove the option "gridSize" in the digitizer workflow O2/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx . Also, we have to check and modify the options initialSpaceChargeDensity (providing file name and histogram name) to calculate the distortion from a histogram and readSpaceCharge (provide a precalculated SpaceCharge object).

Hi @ehellbar,

for the current use in the digitizer we can simply fix the grid size, e.g. to 129 x 129 x 180 (r x z x phi)

Okay, I fixed it to 129x129x180.

The calculations with this grid size are quite fast (~10 sec?)

When I run the external macro to calculate the distortions/corrections for an input histogram it takes ~10sec (8 Threads) for each side.

I replaced the old SpaceCharge class with the new one and it seems that it works in the digitizer. However when I simulate one event and then run the digitizer with
o2-sim-digitizer-workflow -b --distortionType 2 --initialSpaceChargeDensity "/Users/matthias/alice/o2_macros/SpaceCharges/ATO-498/InputSCDensityHistograms_8000events.root,inputSCDensity3D_8000_avg"
it picks up the histogram and calculates the distortions/corrections, but it looks like it performs the calculations multiple times if I dont run the digitizer with --tpc-lanes 1? (In the digitizer it takes around 28 seconds for both sides)

readSpaceCharge (provide a precalculated SpaceCharge object)

So this means the SpaceCharge object is then read from a ROOT file or is it calculated somewhere else and then just passed to the function?

ehellbar · 2020-10-30T17:27:35Z

it picks up the histogram and calculates the distortions/corrections, but it looks like it performs the calculations multiple times if I dont run the digitizer with --tpc-lanes 1? (In the digitizer it takes around 28 seconds for both sides)

Yes, the way this is passed on at the moment, the distortions are calculated for each lane.

So this means the SpaceCharge object is then read from a ROOT file or is it calculated somewhere else and then just passed to the function?

Yes, this is supposed to read a precalculated SpaceCharge object stored in a root file.

ehellbar · 2020-11-02T10:50:54Z

Hello @matthias-kleiner ,
thanks for the updates.

Line 364 in TPCDigitizerSpec.cxx also should be removed:
{"gridSize", VariantType::String, "129,144,129", {"Comma separated list of number of bins in (r,phi,z) for distortion lookup tables (r and z can only be 2**N + 1, N=1,2,3,...)"}},

Can you also adjust the usage in GPU/TPCFastTransformation/macro/generateTPCCorrectionNTuple.C, please? This macro creates the fast correction objects using spline fits of the corrections, to be then applied in the reconstruction.

Cheers,
Ernst

matthias-kleiner · 2020-11-02T11:03:36Z

it picks up the histogram and calculates the distortions/corrections, but it looks like it performs the calculations multiple times if I dont run the digitizer with --tpc-lanes 1? (In the digitizer it takes around 28 seconds for both sides)

Yes, the way this is passed on at the moment, the distortions are calculated for each lane.

So this means the SpaceCharge object is then read from a ROOT file or is it calculated somewhere else and then just passed to the function?

Yes, this is supposed to read a precalculated SpaceCharge object stored in a root file.

I just tried to implement a method to store the SpaceCharge file in a ROOT file, but for some reason it doesnt store the members. I think it might be related to all the static members which are used. But I already implemented a function to set the values (the DataContainer3D member) for the distortions/corrections etc. from a ROOT file. I think it should be easier to store the values for the distortions/corrections in a ROOT file and then set the data members from the ROOT file. Do you think this would also be fine? And if so which values have to be loaded then (only the global distortions and the global corrections)?

ehellbar · 2020-11-02T11:56:53Z

Hi @matthias-kleiner ,
I also had some troubles with the streamers for the old class, and also I just streamed the lookup tables for the global corrections and distortions. So your solution sounds good to me.

matthias-kleiner · 2020-11-02T12:00:34Z

Hello @matthias-kleiner ,
thanks for the updates.

Line 364 in TPCDigitizerSpec.cxx also should be removed:
{"gridSize", VariantType::String, "129,144,129", {"Comma separated list of number of bins in (r,phi,z) for distortion lookup tables (r and z can only be 2**N + 1, N=1,2,3,...)"}},

Can you also adjust the usage in GPU/TPCFastTransformation/macro/generateTPCCorrectionNTuple.C, please? This macro creates the fast correction objects using spline fits of the corrections, to be then applied in the reconstruction.

Cheers,
Ernst

hi @ehellbar,
So if I understand correctly the distortions dont need to be calculated in this macro? So just the global corrections?
And is there some way to test the splines and check if everything is done correctly?
And is there some specific reason the SpaceCharge object is defined as a pointer above the function (o2::tpc::SpaceCharge* sc = 0; ) instead as a simple object in the function o2::tpc::SpaceCharge sc(mNRRows, mNPhiSlices, mNZCols);?

ehellbar · 2020-11-02T13:42:43Z

Hi @matthias-kleiner ,
in the macro generateTPCCorrectionNTuple.C, only the global corrections are actually used and there seems to be no particular reason why o2::tpc::SpaceCharge* sc is a pointer, maybe it is derived from another macro.

But the macro I was actually looking for and which I meant is Detectors/TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C . Here, the AliTPCSpaceCharge3DCalc class is used, so it is not directly affected by the changes, only when we also delete those classes. But it would still be great if you could also change this macro. Also here, the spaceCharge object is a pointer to keep flexibility with grid size and interpolation techniques, but I think in your new implementation this might be redundant and you can make it an object.

There are also two debug functions in createTPCSpaceChargeCorrection.C, one to compare results on the spaceCharge grid and one to compare results on the spline points. In the past, we evaluated the results only for the points for which the corrected point is inside the TPC volume.

matthias-kleiner · 2020-11-02T13:57:15Z

Hi @matthias-kleiner ,
in the macro generateTPCCorrectionNTuple.C, only the global corrections are actually used and there seems to be no particular reason why o2::tpc::SpaceCharge* sc is a pointer, maybe it is derived from another macro.

But the macro I was actually looking for and which I meant is Detectors/TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C . Here, the AliTPCSpaceCharge3DCalc class is used, so it is not directly affected by the changes, only when we also delete those classes. But it would still be great if you could also change this macro. Also here, the spaceCharge object is a pointer to keep flexibility with grid size and interpolation techniques, but I think in your new implementation this might be redundant and you can make it an object.

There are also two debug functions in createTPCSpaceChargeCorrection.C, one to compare results on the spaceCharge grid and one to compare results on the spline points. In the past, we evaluated the results only for the points for which the corrected point is inside the TPC volume.

Hi @ehellbar
thanks for clarification.
I will start looking into the createTPCSpaceChargeCorrection.C macro and delete the old SpaceCharge folder now

matthias-kleiner · 2020-11-03T09:10:07Z

@ehellbar
Is there anything else I should modify before making the pull request?
And do you know why "Check PR with clang-format / build (pull_request_target)" always fails?

ehellbar · 2020-11-03T12:16:49Z

Hi @matthias-kleiner ,
it all looks fine to me.
There are a couple of places where the clang-format failed. You can see the diffs in the "Run clang format" tab after you click on the details.

matthias-kleiner · 2020-11-03T12:25:10Z

Hi @matthias-kleiner ,
it all looks fine to me.
There are a couple of places where the clang-format failed. You can see the diffs in the "Run clang format" tab after you click on the details.

Hi @ehellbar ,
thanks I will look into it and try to fix it.
I should probably also squash all my commits into one commit?

ehellbar · 2020-11-05T09:19:31Z

Hello @matthias-kleiner ,
yes I think it would be good to squash the commit. The conflict with Detectors/TPC/simulation/src/SpaceCharge.cxx should disappear when you apply your commits on a fresh pull.

wiechula

Looks fine. Only two very minor comments.

wiechula · 2020-11-05T10:17:48Z

Detectors/TPC/simulation/src/Digitizer.cxx

Why is SpaceCharge as variable capital?

thanks for noticing it. It should be of course just spaceCharge. I probably changed it by accident

wiechula · 2020-11-05T10:19:21Z

Detectors/TPC/spacecharge/include/TPCSpacecharge/DataContainer3D.h

can this function be constexpr?

yes I changed it

davidrohr · 2020-12-01T16:42:40Z

@matthias-kleiner : Could you apply clang-format :)

davidrohr · 2020-12-01T22:38:08Z

@matthias-kleiner : The macro is now failing in all CI instances, not only MacOS, weird. Perhaps it is related to the problem we had before that some tests were simply not running.
Anyhow, I think the CI right now tests the macro in the interpreted mode. Can it actually run interpreted, or does it need to be compiled? I think if so you have to indicated that in the CMake when you add the macro test.

matthias-kleiner · 2020-12-02T08:35:20Z

@matthias-kleiner : The macro is now failing in all CI instances, not only MacOS, weird. Perhaps it is related to the problem we had before that some tests were simply not running.
Anyhow, I think the CI right now tests the macro in the interpreted mode. Can it actually run interpreted, or does it need to be compiled? I think if so you have to indicated that in the CMake when you add the macro test.

The calculateDistortionsCorrections.C macro needs to be compiled. The other createResidualDistortionObject.C macro should also run in interpreted mode. I will then indicate it with "COMPILE_ONLY"

davidrohr · 2020-12-02T09:31:42Z

@ktf @TimoWilken : Could you check what is the problem with CGAL and GMP. After alisw/alidist#2648 GMP should be a direct build dependency of O2, So how can it be that the GMP include folder does not exist? :

CMake Error in Detectors/TPC/spacecharge/CMakeLists.txt:
  Imported target "CGAL::CGAL" includes non-existent path

    "/jenkins/workspace/DailyBuilds/DailyO2/sw/20156194/1/slc7_x86-64/GMP/v6.0.0-59/include"

I have tried this locally in the slc7-builder container but cannot reproduce it.

TimoWilken · 2020-12-02T09:36:51Z

@davidrohr @ktf That includes something under /jenkins/workspace, but the CI builders are working under /mnt/mesos/sandbox/sandbox. I'd guess there is something wrong with the Jenkins build?

davidrohr · 2020-12-02T14:09:55Z

@matthias-kleiner : Before I was definint the gridsize for the digitizer workflow via --gridSize=65,180,65. This option no longer exists. How do I set the grid size now?

matthias-kleiner · 2020-12-02T14:13:09Z

@matthias-kleiner : Before I was definint the gridsize for the digitizer workflow via --gridSize=65,180,65. This option no longer exists. How do I set the grid size now?

the gridsize is now set to a default of 129,180,129. There is currently no option to set it.

davidrohr · 2020-12-02T14:17:42Z

@matthias-kleiner : I have retried without the gridsize setting, but the digitization output for the TPC looks broken. Below are 2 visualizations of a single pp collision before and after this PR digitized with

o2-sim-digitizer-workflow --TPCtriggered --distortionType 2 --initialSpaceChargeDensity=../InputSCDensityHistograms_8000events.root,inputSCDensity3D_8000_0 --incontext collisioncontext.root --tpc-lanes=1 --onlyDet TPC

using the SCD files provided by Ernst some time ago.

matthias-kleiner · 2020-12-02T14:22:43Z

@matthias-kleiner : I have retried without the gridsize setting, but the digitization output for the TPC looks broken. Below are 2 visualizations of a single pp collision before and after this PR digitized with
o2-sim-digitizer-workflow --TPCtriggered --distortionType 2 --initialSpaceChargeDensity=../InputSCDensityHistograms_8000events.root,inputSCDensity3D_8000_0 --incontext collisioncontext.root --tpc-lanes=1 --onlyDet TPC
using the SCD files provided by Ernst some time ago.

Okay, thanks for doing this check. I will look into this and see if I can figure out what goes wrong

matthias-kleiner · 2020-12-02T15:12:49Z

@davidrohr could you please do the same check again? I just found a bug where I set the z coordinate equal to the y coordinate after the distortion. If this doesnt fix the issue, I could look tomorrow more into it.

davidrohr

@matthias-kleiner : Thx, with this fix it looks OK, and also the reconstruction with the correction looks like before. So if the CI finally passes, we can finally merge.
Could you just quickly comment on the gridSize parameter? I have seen that also in the macro to create the tpcTransform object, the gridSize parameter is gone. So I have just left it out for my tests. Does it not exist anymore, or is the default just ok, or shall I set it somehow else?

matthias-kleiner · 2020-12-02T16:13:00Z

@matthias-kleiner : Thx, with this fix it looks OK, and also the reconstruction with the correction looks like before. So if the CI finally passes, we can finally merge.
Could you just quickly comment on the gridSize parameter? I have seen that also in the macro to create the tpcTransform object, the gridSize parameter is gone. So I have just left it out for my tests. Does it not exist anymore, or is the default just ok, or shall I set it somehow else?

Since the grid size is templated now in the spacecharge class, the easiest solution for now was to just set it to a default value during the simulations and in the macros. Of course the correct default value of the grid size still has to be determined, but gridR=129, gridPhi180, gridZ=129 should already be a good starting point. The grid size is just hardcoded for example in the digitizer using SC = SpaceCharge<double, 129, 129, 180>;. If you want to change the grid size, the only way for now is to change the hardcoded values in Digitizer.h, TPCDigitizerSpec.cxx, and the macros.

Of course it would be better to be able to set the grid size and pass somehow the template parameters of the grid or choose between different presets of grid sizes, but I am not sure if there is an elegant way to implement this in the Digitizer.

sawenzel · 2020-12-02T16:43:47Z

Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceChargeHelpers.h

+  /// analytical potential
+  std::function<DataT(DataT, DataT, DataT)> mPotentialFunc = [& mParA = mParA, &mParB = mParB, &mParC = mParC](const DataT z, const DataT r, const DataT phi) {
+    const DataT zz = std::abs(z);
+    return -mParA * (std::pow((-r + 254.5 + 83.5), 4) - 338.0 * std::pow((-r + 254.5 + 83.5), 3) + 21250.75 * std::pow((-r + 254.5 + 83.5), 2)) * std::cos(mParB * phi) * std::cos(mParB * phi) * std::exp(-1 * mParC * (zz - 125) * (zz - 125));


depending on where or how often this function is used during runtime, I'd like to point out that constructs like std::pow((-r + 254.5 + 83.5), 4) are unfortunately inducing a function call to pow and therefore slow. You can check with https://godbolt.org/ that this is only optimized away with -ffast-math which we are not using.

thanks for the comment :) I didnt knew that about std::pow. But since these functions are only used for "debugging/testing" anyway, it does not really matter how fast the evaluation is.

sawenzel · 2020-12-02T16:46:08Z

Detectors/TPC/spacecharge/src/PoissonSolver.cxx

+  int tnRRow = iOne == 1 ? Nr : Nr / iOne + 1;
+  int tnZColumn = jOne == 1 ? Nz : Nz / jOne + 1;
+
+  std::vector<DataT> coefficient1(Nr);


I am not quite getting the rational why the GRID size is a template argument everywhere. But since it is, there is no need to use dynamic memory in many places. std::array<DataT, Nr> coefficent would be a more sensible choice.

You are right. The grid size could also be just a member in the different classes and then one could use a std::vector everywhere, which would make the code a lot easier. But at the beginning I thought it would make sense to use arrays instead of vectors and to ensure that one uses the correct grid sizes (there are only specific grid sizes which can be used) and that these classes can only be used with the same grid size. I might remove the templates in the future, but this would require some work...

Indeed it makes sense to use an array therestd::array<DataT, Nr> coefficent, but since this vector is passed to several member functions and in a different function in the code the same coefficient array is not depending on the templated grid size std::vector<DataT> coefficient1(tnRRow) it has to be a vector there. I could probably do some additional templating, to use vectors and arrays where appropriate.

sawenzel

It's a nice PR. Please consider my comments for a future polishing step.

davidrohr · 2020-12-02T20:04:20Z

shit, this breaks O2-dataflow defaults, because they disable CGAL.
What do we do. We could remove the required for CGAL, but that will probably be complicated since then we have do disable compiling this code and all that depends on it dynamically depending on whether CGAL is found.
Perhaps I would make it depend on BUILD_SIMULATION:

Add CGAL as simulation dependency in O2SimulationDependencies.cmake.
And then disable the respective libraries when simulation is not present. Might probably be easier since then we don't have to distinguish in the digitizer?

matthias-kleiner · 2020-12-03T12:34:36Z

shit, this breaks O2-dataflow defaults, because they disable CGAL.
What do we do. We could remove the required for CGAL, but that will probably be complicated since then we have do disable compiling this code and all that depends on it dynamically depending on whether CGAL is found.
Perhaps I would make it depend on BUILD_SIMULATION:
* Add CGAL as simulation dependency in O2SimulationDependencies.cmake.

* And then disable the respective libraries when simulation is not present. Might probably be easier since then we don't have to distinguish in the digitizer?

Hi @davidrohr, if the CGAL package causes too much trouble, I could try and see if I can remove it in the space charge class. CGAL is used only in one specific function and I think it might be possible/equally fast to remove CGAL and modify the function where it is used. I can check this and see if its reasonable without CGAL.

davidrohr · 2020-12-03T12:37:39Z

Hi @davidrohr, if the CGAL package causes too much trouble, I could try and see if I can remove it in the space charge class. CGAL is used only in one specific function and I think it might be possible/equally fast to remove CGAL and modify the function where it is used. I can check this and see if its reasonable without CGAL.

OK, can you check that quickly? If it is too much effort, I just discussed with @sawenzel that we could add CGAL as requirement also for o2-dataflow.

davidrohr · 2020-12-03T12:40:47Z

But if that would mean copy&pasting much code, we can also just enforce having CGAL. So perhaps have a quick look what needs to be done and then report back. In the meantime, I am checking what would need to be done to make the CGAL optional.

matthias-kleiner · 2020-12-03T12:47:07Z

But if that would mean copy&pasting much code, we can also just enforce having CGAL. So perhaps have a quick look what needs to be done and then report back. In the meantime, I am checking what would need to be done to make the CGAL optional.

I use CGAL actually only for finding the nearest neighbour of some scattered datapoints and then apply some iterative method to approach some query point. I think i could just skip the nearest neighbour searching and start just at some point in the middle of the grid... There would be actually less code than before, but I have to check how fast it is and if the result will be the same. I can test this later today

davidrohr · 2020-12-03T12:57:02Z

ok, thx, actually we should have thought about that before we implemented the other workarounds...

matthias-kleiner · 2020-12-03T14:59:09Z

ok, thx, actually we should have thought about that before we implemented the other workarounds...

Hi @davidrohr, it looks like I can get rid of the CGAL package and use some slightly different approach in the function where CGAL is used.
The resulting distortions are basically the same and the performance is actually a little bit better.
I am really sorry for all the trouble and work you had with this.

I will then remove all the CGAL stuff from my classes and push it as soon as possible.

davidrohr · 2020-12-03T15:04:28Z

@matthias-kleiner : Ok, great.
And no worries, I guess we all learned something in this endeavor :)

…oximation of nearest neighbour

wiechula requested changes Oct 20, 2020

View reviewed changes

ehellbar reviewed Oct 21, 2020

View reviewed changes

wiechula reviewed Nov 5, 2020

View reviewed changes

matthias-kleiner force-pushed the SpaceChargeDistortionCalculation branch from f37641a to d81a8c4 Compare November 7, 2020 15:44

matthias-kleiner marked this pull request as ready for review November 10, 2020 15:51

matthias-kleiner requested a review from davidrohr as a code owner November 10, 2020 15:51

clang format

5b293b5

changed macro calculateDistortionsCorrections.C to compile only

0a93bf9

fixed z coordinate in distortElectron

1d08bf0

davidrohr previously approved these changes Dec 2, 2020

View reviewed changes

sawenzel reviewed Dec 2, 2020

View reviewed changes

sawenzel previously approved these changes Dec 2, 2020

View reviewed changes

removed CGAL package and replaced nearest neighbour searching by appr…

c60c4da

…oximation of nearest neighbour

matthias-kleiner dismissed stale reviews from sawenzel and davidrohr via c60c4da December 3, 2020 15:19

davidrohr merged commit 55a67d2 into AliceO2Group:dev Dec 3, 2020


		static constexpr unsigned int getID() { return ID; }

		// private:

Adding classes to calculate the space charge distortions and corrections #4618

Adding classes to calculate the space charge distortions and corrections #4618

Uh oh!

Conversation

matthias-kleiner commented Oct 16, 2020

Uh oh!

wiechula left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ehellbar Oct 21, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ehellbar Oct 21, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

ehellbar commented Oct 21, 2020

Uh oh!

matthias-kleiner commented Oct 23, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

matthias-kleiner commented Oct 23, 2020

Uh oh!

ehellbar commented Oct 25, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

matthias-kleiner commented Oct 26, 2020

Uh oh!

matthias-kleiner commented Oct 26, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ehellbar commented Oct 26, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ehellbar commented Oct 26, 2020

Uh oh!

matthias-kleiner commented Oct 30, 2020

Uh oh!

ehellbar commented Oct 30, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

matthias-kleiner commented Oct 30, 2020

Uh oh!

ehellbar commented Oct 30, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ehellbar commented Nov 2, 2020

Uh oh!

matthias-kleiner commented Nov 2, 2020

Uh oh!

ehellbar commented Nov 2, 2020

Uh oh!

matthias-kleiner commented Nov 2, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ehellbar commented Nov 2, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

ehellbar Oct 21, 2020 •

edited

Loading

ehellbar Oct 21, 2020 •

edited

Loading

matthias-kleiner commented Oct 23, 2020 •

edited

Loading

ehellbar commented Oct 25, 2020 •

edited

Loading

matthias-kleiner commented Oct 26, 2020 •

edited

Loading

ehellbar commented Oct 26, 2020 •

edited

Loading

ehellbar commented Oct 30, 2020 •

edited

Loading

ehellbar commented Oct 30, 2020 •

edited

Loading

matthias-kleiner commented Nov 2, 2020 •

edited

Loading

ehellbar commented Nov 2, 2020 •

edited

Loading