Implementation of the algorithm to find the experimental Bragg curves from the AtTracks.

AtData/AtPattern/AtPattern.h

@@ -150,6 +150,20 @@ class AtPattern : public TObject {
    void SetPatternPar(std::vector<double> par) { fPatternPar = std::move(par); }
    void SetChi2(double chi2) { fChi2 = chi2; }
 
+   /**
+    * Calculate the distance in mm along the line from point at parameter t1 to point at parameter t2.
+    * By default returns -99999, which should not be possible in any case.
+    * This function needs to be overriden for any specific AtPattern subclass in order to be useful.
+    */
+   virtual Double_t DistanceAlongPattern(double t1, double t2) const { return -99999; }
+
+   /**
+    * Parameter value at the point passed as input.
+    * By default returns -99999, which should not be possible in any case.
+    * This function needs to be overriden for any specific AtPattern subclass in order to be useful.
+    */
+   virtual double parameterAtPoint(const XYZPoint &point) const { return -99999; }
+
 protected:
    /**
     * Called by other versions of FitPattern. If pointCharge is not empty does charge weighted fit.

AtData/AtPattern/AtPatternLine.cxx

@@ -210,3 +210,12 @@ std::vector<Double_t> AtPatternLine::lineIntersecR(Double_t rMax, Double_t tMin,
    }
    return result;
 }
+
+Double_t AtPatternLine::DistanceAlongPattern(double t1, double t2) const
+{
+   XYZPoint point1 = GetPointAt(t1);
+   XYZPoint point2 = GetPointAt(t2);
+   Double_t distanceSquared = std::pow(point1.X() - point2.X(), 2) + std::pow(point1.Y() - point2.Y(), 2) +
+                              std::pow(point1.Z() - point2.Z(), 2);
+   return TMath::Sqrt(distanceSquared);
+}

AtData/AtPattern/AtPatternLine.h

@@ -43,11 +43,13 @@ class AtPatternLine : public AtPattern {
 
    TEveLine *GetEveLine(Double_t rMax = 250) const;
 
+   virtual Double_t DistanceAlongPattern(double t1, double t2) const override;
+   double parameterAtPoint(const XYZPoint &point) const override;
+
 protected:
    std::vector<Double_t> lineIntersecR(Double_t rMax, Double_t tMin, Double_t tMax) const;
 
    virtual void FitPattern(const std::vector<XYZPoint> &points, const std::vector<double> &charge) override;
-   double parameterAtPoint(const XYZPoint &point) const;
    ClassDefOverride(AtPatternLine, 1)
 };
 } // namespace AtPatterns

AtData/AtTrack.cxx

@@ -28,7 +28,7 @@ AtTrack &AtTrack::operator=(AtTrack obj)
 AtTrack::AtTrack(const AtTrack &o)
    : fTrackID(o.fTrackID), fIsMerged(o.fIsMerged), fVertexToZDist(o.fVertexToZDist), fGeoThetaAngle(o.fGeoThetaAngle),
      fGeoPhiAngle(o.fGeoPhiAngle), fGeoRadius(o.fGeoRadius), fGeoCenter(o.fGeoCenter),
-     fHitClusterArray(o.fHitClusterArray)
+     fHitClusterArray(o.fHitClusterArray), fBraggCurveValues(o.fBraggCurveValues), fBraggCurve(o.fBraggCurve)
 {
    fPattern = (o.fPattern != nullptr) ? o.fPattern->Clone() : nullptr;
    for (auto &hit : o.fHitArray)

AtData/AtTrack.h

@@ -24,6 +24,16 @@ class TMemberInspector;
 
 class AtTrack : public TObject {
 
+public:
+   struct BraggCurve {
+      std::vector<Double_t> IntegratedELossValues;
+      std::vector<Double_t> RangeValues;
+      std::vector<Double_t> ELossErrors;
+      Int_t nBins{0};
+      Double_t binSize{0};
+      Int_t smoothingSteps{0};
+   };
+
 protected:
    using XYZPoint = ROOT::Math::XYZPoint;
    using HitPtr = std::unique_ptr<AtHit>;
@@ -44,6 +54,12 @@ class AtTrack : public TObject {
    std::pair<Double_t, Double_t> fGeoCenter; //< Center of the spiral track
    std::vector<AtHitCluster> fHitClusterArray; //< Clusterized hits container. Can also be stored in fHitArray
 
+   // Obtained in AtPatternModification
+   // Vector of pair of values for the Bragg curve of the track (archLength[mm], eLoss[a.u.]).
+   std::vector<std::pair<Double_t, Double_t>> fBraggCurveValues;
+   // Container for the Bragg curve integrated over the binning.
+   BraggCurve fBraggCurve;
+
 public:
    AtTrack() = default;
    AtTrack(const AtTrack &obj);
@@ -72,6 +88,9 @@ class AtTrack : public TObject {
       swap(a.fGeoPhiAngle, b.fGeoPhiAngle);
       swap(a.fGeoRadius, b.fGeoRadius);
       swap(a.fGeoCenter, b.fGeoCenter);
+
+      swap(a.fBraggCurveValues, b.fBraggCurveValues);
+      swap(a.fBraggCurve, b.fBraggCurve);
    };
 
    // Getters
@@ -89,6 +108,9 @@ class AtTrack : public TObject {
    std::pair<Double_t, Double_t> GetGeoCenter() const { return fGeoCenter; }
    std::vector<AtHitCluster> *GetHitClusterArray() { return &fHitClusterArray; }
 
+   std::vector<std::pair<Double_t, Double_t>> GetBraggCurveValues() { return fBraggCurveValues; }
+   const BraggCurve GetBraggCurve() { return fBraggCurve; }
+
    Bool_t GetIsMerged() const { return fIsMerged; }
    Double_t GetVertexToZDist() const { return fVertexToZDist; }
 
@@ -104,6 +126,9 @@ class AtTrack : public TObject {
    void SetGeoCenter(std::pair<Double_t, Double_t> center) { fGeoCenter = center; }
    void AddClusterHit(std::shared_ptr<AtHitCluster> hitCluster);
 
+   void AddBraggCurvePair(Double_t range, Double_t eLoss) { fBraggCurveValues.push_back(std::make_pair(range, eLoss)); }
+   void SetBraggCurve(BraggCurve braggCurve) { fBraggCurve = braggCurve; }
+
    void SetIsMerged(bool val) { fIsMerged = val; }
    void SetVertexToZDist(Double_t val) { fVertexToZDist = val; }
 
@@ -139,7 +164,7 @@ class AtTrack : public TObject {
       return o;
    }
 
-   ClassDef(AtTrack, 3);
+   ClassDef(AtTrack, 4);
 };
 
 #endif

AtEventDisplay/AtEventDisplayLinkDef.h

@@ -20,6 +20,7 @@
 #pragma link C++ class AtTabMacro + ;
 #pragma link C++ class AtTabEnergyLoss + ;
 #pragma link C++ class AtTabFF + ;
+#pragma link C++ class AtTabBraggCurve + ;
 
 #pragma link C++ class AtTabInfoBase - !;
 #pragma link C++ class AtTabInfo - !;

AtEventDisplay/AtTabs/AtTabBraggCurve.cxx

@@ -0,0 +1,169 @@
+#include "AtTabBraggCurve.h"
+
+#include "AtPattern.h" // for AtPattern
+#include "AtPatternEvent.h"
+#include "AtTabInfo.h" // for AtTabInfoFairRoot, AtTabInfo
+#include "AtTrack.h"   // for AtTrack
+
+#include <FairLogger.h> // for LOG
+
+#include <TAttMarker.h> // for TAttMarker
+#include <TCanvas.h>
+#include <TEveBrowser.h>
+#include <TEveElement.h>      // for TEveElement
+#include <TEveEventManager.h> // for TEveEventManager
+#include <TEveGeoNode.h>
+#include <TEveManager.h>  // for TEveManager, gEve
+#include <TEvePointSet.h> // for TEvePointSet
+#include <TEveViewer.h>
+#include <TEveWindow.h>
+#include <TGLViewer.h>
+#include <TGTab.h>
+#include <TGeoManager.h>
+#include <TRootEmbeddedCanvas.h>
+
+#include <array>   // for array
+#include <utility> // for move
+namespace DataHandling {
+class AtSubject;
+}
+
+ClassImp(AtTabBraggCurve);
+
+AtTabBraggCurve::AtTabBraggCurve() : AtTabMain() {}
+
+AtTabBraggCurve::~AtTabBraggCurve()
+{
+   fPadNum->Detach(this);
+   fEventBranch->Detach(this);
+   fRawEventBranch->Detach(this);
+   fPatternEventBranch->Detach(this);
+   fEntry->Detach(this);
+}
+
+void AtTabBraggCurve::Update(DataHandling::AtSubject *sub)
+{
+   // If we should update the stuff that depends on the AtEvent
+   if (sub == fEventBranch || sub == fEntry) {
+      UpdateEventElements();
+   }
+   if (sub == fPatternEventBranch || sub == fEntry) {
+      UpdatePatternEventElements();
+   }
+
+   // If we should update the 3D display
+   if (sub == fEventBranch || sub == fPatternEventBranch || sub == fEntry) {
+      gEve->Redraw3D(false); // false -> don't reset camera
+   }
+}
+
+void AtTabBraggCurve::MakeTab(TEveWindowSlot *slot)
+{
+   TEveWindowPack *pack = nullptr;
+
+   // 3D
+   pack = slot->MakePack();
+   pack->SetElementName("BraggCurve");
+   pack->SetHorizontal();
+   pack->SetShowTitleBar(kFALSE);
+
+   pack->NewSlot()->MakeCurrent();
+   TEveViewer *view3D = gEve->SpawnNewViewer("3D View", "");
+   view3D->AddScene(gEve->GetGlobalScene());
+   view3D->AddScene(gEve->GetEventScene());
+
+   slot = pack->NewSlot();
+   TEveWindowPack *pack2 = slot->MakePack();
+   pack2->SetShowTitleBar(kFALSE);
+   pack2->SetVertical();
+   slot = pack2->NewSlot();
+   slot->StartEmbedding();
+   // fCvsDeDx = new TCanvas("dEdx Bragg curve Canvas");
+   // fCvsDeDx->ToggleEditor();
+   slot->StopEmbedding();
+
+   slot = pack2->NewSlotWithWeight(1.5);
+   auto *ecvs = new TRootEmbeddedCanvas();
+   TEveWindowFrame *frame = slot->MakeFrame(ecvs);
+   frame->SetElementName("Bragg curve Canvas");
+   pack->GetEveFrame()->SetShowTitleBar(kFALSE);
+   fCvsELossVRange = ecvs->GetCanvas();
+   // fCvsELossVRange->AddExec("ex", "AtTab3DBraggCurve::NextTrack()");
+
+   fCvsELossVRange->ToggleEventStatus();
+   DrawHistELossVRange();
+
+   if (gGeoManager) {
+      TGeoNode *geoNode = gGeoManager->GetTopNode();
+      Int_t option = 1;
+      Int_t level = 3;
+      Int_t nNodes = 10000;
+      auto *topNode = new TEveGeoTopNode(gGeoManager, geoNode, option, level, nNodes);
+      gEve->AddGlobalElement(topNode);
+
+      Int_t transparency = 80;
+      gGeoManager->GetVolume("drift_volume")->SetTransparency(transparency);
+      gEve->FullRedraw3D(kTRUE);
+   }
+
+   gEve->GetBrowser()->GetTabRight()->SetTab(1);
+
+   gEve->Redraw3D(true, true);
+
+   TGLViewer *dfViewer = gEve->GetDefaultGLViewer(); // Is this doing anything?
+   dfViewer->CurrentCamera().RotateRad(-.7, 0.5);
+   dfViewer->DoDraw();
+   UpdateRenderState();
+}
+
+void AtTabBraggCurve::DrawHistELossVRange()
+{
+   AtTrack::BraggCurve braggCurve;
+   braggCurve.nBins = 1000;
+   braggCurve.binSize = 1;
+   DrawHistELossVRange(braggCurve);
+}
+
+void AtTabBraggCurve::DrawHistELossVRange(AtTrack::BraggCurve braggCurve)
+{
+   if (fHistELossVRange != nullptr)
+      fCvsELossVRange->GetListOfPrimitives()->Remove(fHistELossVRange);
+
+   int nBins = braggCurve.nBins;
+   double binSize = braggCurve.binSize;
+
+   fHistELossVRange = new TH1F("Charge vs Range", "Charge vs Range", nBins, 0, nBins * binSize);
+   fHistELossVRange->SetDirectory(0);
+   fCvsELossVRange->cd();
+   fHistELossVRange->Draw();
+   fHistELossVRange->GetXaxis()->SetTitle("Range [mm]");
+   fHistELossVRange->GetYaxis()->SetTitle("Charge [ADC]");
+
+   for (int i = 0; i < braggCurve.IntegratedELossValues.size(); i++) {
+      fHistELossVRange->SetBinContent(i + 1, braggCurve.IntegratedELossValues[i]);
+      fHistELossVRange->SetBinError(i + 1, braggCurve.ELossErrors[i]);
+   }
+
+   fCvsELossVRange->Modified();
+   fCvsELossVRange->Update();
+}
+
+void AtTabBraggCurve::UpdatePatternEventElements()
+{
+   AtTabMain::UpdatePatternEventElements();
+
+   auto fPatternEvent = GetFairRootInfo<AtPatternEvent>();
+   if (fPatternEvent == nullptr) {
+      LOG(debug) << "Cannot update AtPatternEvent elements: no event available";
+      return;
+   }
+
+   auto &tracks = fPatternEvent->GetTrackCand();
+   if (tracks.size()) {
+      fTrackIdx = 0;
+      DrawHistELossVRange(tracks[fTrackIdx].GetBraggCurve());
+   } else {
+      fTrackIdx = -1;
+      DrawHistELossVRange();
+   }
+}

AtEventDisplay/AtTabs/AtTabBraggCurve.h

@@ -0,0 +1,49 @@
+#ifndef ATTABBRAGGCURVE_H
+#define ATTABBRAGGCURVE_H
+
+#include "AtTabMain.h" // for AtTabMain::TEvePointSetPtr, AtTa...
+#include "AtTrack.h"
+#include "AtViewerManagerSubject.h" // for AtBranch
+
+#include <Rtypes.h>           // for THashConsistencyHolder, ClassDef...
+#include <TEveEventManager.h> // for TEveEventManager
+#include <TEvePointSet.h>     // for TEvePointSet
+#include <TH1F.h>
+
+#include <array>        // for array
+#include <memory>       // for make_unique
+class TBuffer;          // lines 21-21
+class TClass;           // lines 23-23
+class TMemberInspector; // lines 27-27
+namespace DataHandling {
+class AtSubject;
+}
+
+/**
+ * @brief Tab for visualizing the Bragg curves of an AtTrack
+ */
+class AtTabBraggCurve : public AtTabMain {
+protected:
+   TCanvas *fCvsELossVRange{nullptr};
+   TH1F *fHistELossVRange{nullptr};
+
+   int fTrackIdx{-1};
+
+public:
+   AtTabBraggCurve();
+   ~AtTabBraggCurve();
+
+   virtual void Update(DataHandling::AtSubject *sub) override;
+
+protected:
+   virtual void MakeTab(TEveWindowSlot *slot) override;
+
+   void DrawHistELossVRange();
+   void DrawHistELossVRange(AtTrack::BraggCurve braggCurve);
+
+private:
+   virtual void UpdatePatternEventElements() override;
+
+   ClassDefOverride(AtTabBraggCurve, 1);
+};
+#endif

AtEventDisplay/AtTabs/AtTabMain.h

@@ -95,7 +95,7 @@ class AtTabMain : public AtTabBase, public DataHandling::AtObserver {
    void SetPointsFromHits(TEvePointSet &hitSet, const std::vector<AtHit *> &hits);
    void SetPointsFromTrack(TEvePointSet &hitSet, const AtTrack &track);
 
-private:
+protected:
    // Functions to draw the initial canvases
    void DrawPadPlane();
    void DrawPadWave();
@@ -105,7 +105,7 @@ class AtTabMain : public AtTabBase, public DataHandling::AtObserver {
    // Update hit sets
    void UpdatePadPlane();
    void UpdateEventElements();
-   void UpdatePatternEventElements();
+   virtual void UpdatePatternEventElements();
 
    void ExpandNumPatterns(int num);
 

AtEventDisplay/CMakeLists.txt

@@ -26,6 +26,7 @@ AtTabs/AtTabPad.cxx
 AtTabs/AtTabInfoTree.cxx
 AtTabs/AtTabEnergyLoss.cxx
 AtTabs/AtTabFF.cxx
+AtTabs/AtTabBraggCurve.cxx
 
 AtSidebar/AtEventSidebar.cxx
 AtSidebar/AtSidebarFrames.cxx