Implemented the 3D Bragg curve fitter using CATima by RealAurio · Pull Request #250 · ATTPC/ATTPCROOTv2

RealAurio · 2025-06-10T08:06:10Z

Been working on this feature. Finally managed to get something that seems to work. There is still a lot to fix, and probably a refactoring will be done once I have everything more clear. However, I wanted to create this pull request so that we can check how does this feature performs with the simulations that are being prepared by others in their respective forks.

CATima repo: https://github.com/hrosiak/catima

Notes:

In order to get this to work, I needed to compile CATima using the DGSL integration option (cmake -DCMAKE_INSTALL_PREFIX=~/fair_install/catimaInstall -DGSL_INTEGRATION=ON ..).
I also needed to comment 2 lines in the file "src/storage.cpp" file from CATima (in deleter function InterpolatorGSL::~InterpolatorGSL(), just comment both lines). Maybe it was just me, but it was giving me a segmentation fault.
Once you have installed CATima, export the environment variable CATIMA and set it to the path where you have installed CATima (it should have 3 folders, namely: include, lib and share). In my case: export CATIMA=~/fair_install/catimaInstall
Now simply compile ATTPCROOTv2. It should find CATima and therefore compile the At3DBraggFitter class.
You can test that it works by running the simulation in "macro/Simulation/ATTPC/12Be_dp".
I have removed GenFit from the AtFitter base class, since I checked it and it really doesn't use anything from GenFit. I initially wanted At3DBraggFitter to derive from AtFitter, so that is why I made that change. I hope it doesn't break anything.
Also, there are a couple of changes regarding the S800, since I was trying to test this with the data from the e23031 experiment. Sorry for mixing things.

…he e23031 experiment, which eventually will be used to test the PID with real data.

…n, and as a means of PID for straight-tracks.

anthoak13

The biggest issue that needs fixing is the fact it can't compile without catima!

Other high-level notes:

New tab re-implements a lot of AtTabMain. Should either inherit, or just be a normal tab.
Catima is tightly connected to the bragg class. Should be integrated as an AtELossModel and that is what is exposed to the Brag classes so we have a consistent handling of energy loss models throughout the code, and to make catima integration useful for things like simulation and other fitters (UKF, MCFitter, etc).
I don't understand the purpose of the handler class, it seems necessary?
This class should NOT inherit from AtFitter as written. It breaks those conventions
We should think about fit results and merge the different fit result data classes.
Data flow (copying tracks, etc. seems fishy. We should thing about that more (related to prev point)

anthoak13 · 2025-06-10T21:45:52Z

AtData/At3DBraggFitResult.h

+#include <utility>
+#include <vector>
+
+class At3DBraggFitResult : public TObject {


This class is similar in scope to the AtMCResult class that exists. I think we should talk about how to merge these two - maybe this derives from AtMCResult and that gets refactored to AtFitResult?

This and the MC fitting are basically doing the same thing - only one more than once per species.

I didn't know about that class. I will take a look into it.

anthoak13 · 2025-06-10T21:46:45Z

AtData/AtTrack.h

+   // Stuff related to the 3DBragg curve. May be moved to a child class in the future.
+   // Vector of pair of values for the 3DBragg curve of the track (archLength[mm], eLoss[a.u.]).
+   std::vector<std::pair<Double_t, Double_t>> f3DBraggCurveValues;
+   std::unique_ptr<At3DBraggFitResult> fFitResult{};


I don't like this being in the AtTrack class given it's specific to one of the three fitting methods we posses.

As a first thought, I would have the fit result as a separate branch on the tree. You can just save the FitResult class and store the track ID there for linking.

Yeah, that makes much more sense. I will probably start with this once I start fixing things.

anthoak13 · 2025-06-10T21:49:48Z

AtData/AtTrack.h


+   // Stuff related to the 3DBragg curve. May be moved to a child class in the future.
+   // Vector of pair of values for the 3DBragg curve of the track (archLength[mm], eLoss[a.u.]).
+   std::vector<std::pair<Double_t, Double_t>> f3DBraggCurveValues;


Is this just a copy of the fit result?

f3DBraggCurveValues is not the fit result, but rather the experimental values of (ELoss, range) to be fitted. The fit result is in line 51 (fFitResult).

As I initially implemented this, first the At3DBraggCurveHandler class takes the AtTrack and just obtains these (ELoss, range) by projecting the trace of each AtHit into the RANSAC AtPatternLine. Then, after these pairs are obtained, the fitter class is used in order to fit them to a Bragg curve.

My thought in separating this process in these 2 steps is that one can first only use the At3DBraggCurveHandler to check in what order the tail of the Bragg curve is in ADC, so the amplitude factor of the fitter can be tuned, and then repeat the analysis with the fitter on.

Also, since this will eventually be refactored so that it is stored in its own branch, both the (ELoss, range) pairs and the fitting result will be stored in a new AtBraggCurve class object in the AtData folder, maybe?

anthoak13 · 2025-06-10T21:54:01Z

AtEventDisplay/AtTabs/AtTab3DBraggCurve.h

+   Double_t fDZ{1};
+   Double_t fDS{1};
+
+   std::unique_ptr<AtFITTER::At3DBraggFitter> fFitter{nullptr};


Why does the tab have ownership over the fitter object?

It's just so that the best fitting dEdx curves for each candidate particle can be calculated and plotted on the top-right window of this tab (I attach a screenshot).

Initially, I was calculating these curves in the At3DBraggFitter class, and then storing them in the At3DBraggFitResult class. However, this made the output files too heavy due to the amount of points needed to have good resolution of the peak, and the analysis was around 4 times slower. So eventually, I just decided to calculate them in the tab with the stored fitted kinetic energy. This makes that panel quite slow and laggy tho, so it probably needs to be improved.

anthoak13 · 2025-06-10T21:55:33Z

AtEventDisplay/AtTabs/AtTab3DBraggCurve.h

+
+   void SetDZ(Double_t DZ) { fDZ = DZ; }
+
+   void Set3DBraggCurveFitter(std::unique_ptr<AtFITTER::At3DBraggFitter> fitter) { fFitter = std::move(fitter); }


If the fitter is required for the tab to function, it would make more sense to pass it in the constructor

Well, the tab can function without the fitter. If no fitter is provided, then the dEdx curves on the top panel are not plotted. Since right now it's quite laggy, I thought it was a good idea to leave the option to not provide a fitter, which prevents the dEdx curves from being calculated and plotted.

anthoak13 · 2025-06-10T23:03:50Z

AtReconstruction/AtFitter/At3DBraggFitter.cxx

+   return dEdxVector;
+}
+
+std::vector<std::pair<double, double>> AtFITTER::At3DBraggFitter::GetELossModel(Double_t kineticEnergy)


What is this? How is it different than the previous method?

My brain is staring to get a little tired, so I might have more questions lol

Yeah, I should have commented.

The previous function computes the dEdx Vs range curve, with "good resolution". However, this function does essentially the same thing, but integrates dEdx in each bin of the histogram. Since our experimental data is the charge deposited in each bin, we need to compare that with the energy deposited in that bin, which is the integral of dEdx. In the lower panel of the tab, the black curve is the ELoss (integral of dEdx) vs range values obtained with this method, thus the low resolution.

Initially, I only had this method, and then decided to compute the dEdx curves so they could be shown, which is why both methods seem to be doing the same thing. However, I should refactor this so that the GetELossModel calls the GetDeDx, and then integrates de dE/dx bin by bin.

anthoak13 · 2025-06-10T23:07:34Z

AtReconstruction/AtPatternRecognition/AtSampleConsensus.cxx

+   if (auto patternLine = dynamic_cast<AtPatterns::AtPatternLine *>(pattern)) {
+      track.SetGeoTheta(TMath::ATan2(
+         std::abs(patternLine->GetDirection().Z()),
+         TMath::Sqrt(std::pow(patternLine->GetDirection().X(), 2) + std::pow(patternLine->GetDirection().Y(), 2))));
+      track.SetGeoPhi(TMath::ATan2(patternLine->GetDirection().Y(), patternLine->GetDirection().X()));
+   }
+


I think we require c++17 so fine, but we should probably double check our min version isn't 14

I think I don't understand this

anthoak13 · 2025-06-10T23:09:23Z

AtReconstruction/AtPulseAnalyzer/AtPSA.h

+   std::uniform_real_distribution<double> fUniform{-0.5, 0.5};
+   std::mt19937_64 fRNG;


I don't like these here. They definitely can't stay as is. They are not inherently thread safe (and PSA methods are called in parallel by parts of the code). Should probably be changed to thread local, and moved to the specific PSA you want them for.

Yeah, I should move them to the AtPSAHitPerTBInRegion. I think initially I had a problem with this, so that's why I moved them here, but I will try to move them back.

anthoak13 · 2025-06-10T23:12:59Z

AtReconstruction/AtPulseAnalyzer/AtPSAHitPerTBInRegion.cxx

+void AtPSAHitPerTBInRegion::SetTBLimits(std::pair<Int_t, Int_t> limits)
+{
+   if (limits.first >= limits.second) {
+      std::cout << " Warning AtPSAHitPerTBInRegion::SetTBLimits -  Wrong Time Bucket limits. Setting default limits "
+                   "(0,512) ... "
+                << "\n";
+      fIniTB = 0;
+      fEndTB = 512;
+
+   } else {
+      if (limits.first < 0)
+         fIniTB = 0;
+      else
+         fIniTB = limits.first;
+
+      if (limits.second > fNumTbs)
+         fEndTB = fNumTbs;
+      else
+         fEndTB = limits.second;
+   }
+}


This looks like just a better version of the base method. As such, it should just replace the base method - not be overwritten

It may be a very noob question but, what's the difference?

Ahhhhhhh, I just understood. I am currently redoing everything in this pull request, and I am at this point. I will replace the method in AtPSAHitPerTB directly.
In fact, this whole class only differs on AtPSAHitPerTB in one detail. The base class, after all subHits have been generated, it sets their trace integral to the trace integral of the entire hit, while this subclass sets the trace integral of each subHit to the ADC content of the bin corresponding to it.
I think I will not create a new class, I will simply edit the base one, adding a bool member which controls the setting of the trace integral at the end (set to true by default so nothing that may be using this breaks).

anthoak13 · 2025-06-10T23:13:51Z

AtReconstruction/CMakeLists.txt

  AtPatternRecognition/triplclust/src/directedgraph.cxx
  AtPatternRecognition/triplclust/src/postprocess.cxx
  AtPatternRecognition/AtTrackFinderTC.cxx
+  AtPatternRecognition/At3DBraggCurveHandler.cxx


I think something here depends on catima

At3DBraggCurveHandler doesn't depend on CATima. However, since I included At3DBraggFitter in the At3DBraggCurveTask and in the AtTab3DBraggCurve without thinking, the dependence on CATima is in those classes by accident.

… matching.

anthoak13 · 2025-08-05T18:30:01Z

Replaced by #257

RealAurio added 2 commits June 10, 2025 16:16

Implemented the E1_dE from the S800 and added the macros related to t…

0adc392

…he e23031 experiment, which eventually will be used to test the PID with real data.

Implemented the fitting of the 3D Bragg curve for energy recontructio…

071603d

…n, and as a means of PID for straight-tracks.

anthoak13 requested changes Jun 10, 2025

View reviewed changes

RealAurio added 2 commits June 12, 2025 15:53

Divided the charge of the big pads by 4 to try to compensate the gain…

fb9432c

… matching.

Why? Just why?

305c666

anthoak13 closed this Aug 5, 2025


		void SetDZ(Double_t DZ) { fDZ = DZ; }

		void Set3DBraggCurveFitter(std::unique_ptr<AtFITTER::At3DBraggFitter> fitter) { fFitter = std::move(fitter); }

		std::uniform_real_distribution<double> fUniform{-0.5, 0.5};
		std::mt19937_64 fRNG;

Comments

Conversation

RealAurio commented Jun 10, 2025

Uh oh!

anthoak13 left a comment • 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!

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!

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!

Choose a reason for hiding this comment

Uh oh!

anthoak13 commented Aug 5, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

anthoak13 left a comment •

edited

Loading