Adding observables that measure the sigma of the distance of the hits to the line

inc/TRestTrackEvent.h

@@ -88,6 +88,7 @@ class TRestTrackEvent : public TRestEvent {
     Double_t GetEnergy(TString option = "");
 
     Int_t GetLevel(Int_t tck);
+    Int_t GetLevelById(Int_t tckId);
     void SetLevels();
     inline Int_t GetLevels() const { return fLevels; }
 

inc/TRestTrackLineAnalysisProcess.h

@@ -57,6 +57,8 @@ class TRestTrackLineAnalysisProcess : public TRestEventProcess {
 
     const char* GetProcessName() const override { return "trackLineAna"; }
 
+    TVector3 GetSigmaToLine(TRestTrack* track, TRestTrack* line, bool ponderateByEnergy = true);
+
     // Constructor
     TRestTrackLineAnalysisProcess();
     // Destructor

src/TRestTrackEvent.cxx

@@ -262,6 +262,22 @@ Int_t TRestTrackEvent::GetLevel(Int_t tck) {
     return lvl;
 }
 
+Int_t TRestTrackEvent::GetLevelById(Int_t tckId) {
+    Int_t lvl = 1;
+    auto track = GetTrackById(tckId);
+    if (track == nullptr) {
+        RESTWarning << "Track with ID " << tckId << " not found" << RESTendl;
+        return -1;
+    }
+    Int_t parentTrackId = track->GetParentID();
+
+    while (parentTrackId > 0) {
+        lvl++;
+        parentTrackId = GetTrackById(parentTrackId)->GetParentID();
+    }
+    return lvl;
+}
+
 Bool_t TRestTrackEvent::isTopLevel(Int_t tck) {
     if (GetLevels() == GetLevel(tck)) return true;
     return false;

src/TRestTrackLineAnalysisProcess.cxx

@@ -214,6 +214,54 @@ TRestEvent* TRestTrackLineAnalysisProcess::ProcessEvent(TRestEvent* inputEvent)
     fOutTrackEvent->AddTrack(tckY);
 
     fOutTrackEvent->SetLevels();
+
+    // Get the track which contains the hits from where the line track was obtained
+    // for the XZ projection
+    auto originalTrackX = fOutTrackEvent->GetTrackById(tckX->GetParentID());
+    while (fOutTrackEvent->GetLevelById(originalTrackX->GetTrackID()) > 1) {
+        originalTrackX = fOutTrackEvent->GetTrackById(originalTrackX->GetParentID());
+    }
+    // for the YZ projection
+    auto originalTrackY = fOutTrackEvent->GetTrackById(tckY->GetParentID());
+    while (fOutTrackEvent->GetLevelById(originalTrackY->GetTrackID()) > 1) {
+        originalTrackY = fOutTrackEvent->GetTrackById(originalTrackY->GetParentID());
+    }
+    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
+        fOutTrackEvent->PrintEvent();
+    }
+    RESTDebug << "Original track X ID: " << originalTrackX->GetTrackID()
+              << "; line track ID: " << tckX->GetTrackID() << RESTendl;
+    RESTDebug << "Original track Y ID: " << originalTrackY->GetTrackID()
+              << "; line track ID: " << tckY->GetTrackID() << RESTendl;
+
+    auto sigmaXZ = GetSigmaToLine(originalTrackX, tckX, true);
+    auto sigmaYZ = GetSigmaToLine(originalTrackY, tckY, true);
+    auto meanSigmaZ = TMath::Sqrt((sigmaXZ.Z() * sigmaXZ.Z() + sigmaYZ.Z() * sigmaYZ.Z()) * 0.5);
+    double totalSigma =
+        TMath::Sqrt(sigmaXZ.X() * sigmaXZ.X() + sigmaYZ.Y() * sigmaYZ.Y() + meanSigmaZ * meanSigmaZ);
+
+    SetObservableValue("energySigmaX", sigmaXZ.X());
+    SetObservableValue("energySigmaY", sigmaYZ.Y());
+    SetObservableValue("energySigmaZ", meanSigmaZ);
+    SetObservableValue("energySigmaXZ", sigmaXZ.Mag());
+    SetObservableValue("energySigmaYZ", sigmaYZ.Mag());
+    SetObservableValue("energySigmaZ_XZ", sigmaXZ.Z());
+    SetObservableValue("energySigmaZ_YZ", sigmaYZ.Z());
+    SetObservableValue("energySigma", totalSigma);
+
+    sigmaXZ = GetSigmaToLine(originalTrackX, tckX, false);
+    sigmaYZ = GetSigmaToLine(originalTrackY, tckY, false);
+    meanSigmaZ = TMath::Sqrt((sigmaXZ.Z() * sigmaXZ.Z() + sigmaYZ.Z() * sigmaYZ.Z()) * 0.5);
+    totalSigma = TMath::Sqrt(sigmaXZ.X() * sigmaXZ.X() + sigmaYZ.Y() * sigmaYZ.Y() + meanSigmaZ * meanSigmaZ);
+    SetObservableValue("sigmaX", sigmaXZ.X());
+    SetObservableValue("sigmaY", sigmaYZ.Y());
+    SetObservableValue("sigmaZ", meanSigmaZ);
+    SetObservableValue("sigmaXZ", sigmaXZ.Mag());
+    SetObservableValue("sigmaYZ", sigmaYZ.Mag());
+    SetObservableValue("sigmaZ_XZ", sigmaXZ.Z());
+    SetObservableValue("sigmaZ_YZ", sigmaYZ.Z());
+    SetObservableValue("sigma", totalSigma);
+
     return fOutTrackEvent;
 }
 
@@ -222,3 +270,65 @@ TRestEvent* TRestTrackLineAnalysisProcess::ProcessEvent(TRestEvent* inputEvent)
 /// processed. Nothing to do...
 ///
 void TRestTrackLineAnalysisProcess::EndProcess() {}
+
+///////////////////////////////////////////////
+/// \brief Function to calculate the sigma of the distances of the hits to the line
+/// defined by the line track. The sigma is calculated as the standard deviation of the
+/// distances of the hits to the line. The line is defined by the two closest nodes to the
+/// hit. The sigma can be ponderated by the energy of the hit or not. This function works
+/// for any type of track (XZ, YZ or XYZ).
+/// TODO: should this function be in TRestVolumeHits ?
+/// \param track The track containing the hits to calculate the sigma
+/// \param line The nodes that define the line.
+/// \param ponderateByEnergy If true, the sigma is ponderated by the energy of the hit
+/// \return The sigma of the distances of the hits to the line
+///
+TVector3 TRestTrackLineAnalysisProcess::GetSigmaToLine(TRestTrack* track, TRestTrack* line,
+                                                       bool ponderateByEnergy) {
+    TRestVolumeHits* lineVolHits = line->GetVolumeHits();
+    TVector3 sigma2ToLine = TVector3(0, 0, 0);
+    for (size_t i = 0; i < track->GetVolumeHits()->GetNumberOfHits(); i++) {
+        TVector3 hitPos = track->GetVolumeHits()->GetPosition(i);
+        auto hitEnergy = track->GetVolumeHits()->GetEnergy(i);
+        size_t closestNodeIndex = lineVolHits->GetClosestHit(hitPos);
+
+        // find second closes node
+        size_t secondClosestNodeIndex = closestNodeIndex - 1;
+        if (secondClosestNodeIndex < 0) {
+            secondClosestNodeIndex = closestNodeIndex + 1;
+        } else if (closestNodeIndex == lineVolHits->GetNumberOfHits() - 1) {
+            secondClosestNodeIndex = closestNodeIndex - 1;
+        } else {
+            auto dist1 = (lineVolHits->GetPosition(secondClosestNodeIndex) - hitPos).Mag();
+            auto dist2 = (lineVolHits->GetPosition(closestNodeIndex + 1) - hitPos).Mag();
+            if (dist1 > dist2) {
+                secondClosestNodeIndex = closestNodeIndex + 1;
+            }
+        }
+        TVector3 lineDirector =
+            lineVolHits->GetPosition(closestNodeIndex) - lineVolHits->GetPosition(secondClosestNodeIndex);
+        TVector3 hitToClosestLinePoint = hitPos - lineVolHits->GetPosition(closestNodeIndex);
+        TVector3 hitToLine =
+            hitToClosestLinePoint - (hitToClosestLinePoint.Dot(lineDirector.Unit())) * lineDirector.Unit();
+
+        double toAddX = hitToLine.X() * hitToLine.X();
+        double toAddY = hitToLine.Y() * hitToLine.Y();
+        double toAddZ = hitToLine.Z() * hitToLine.Z();
+        if (ponderateByEnergy) {
+            toAddX *= hitEnergy;
+            toAddY *= hitEnergy;
+            toAddZ *= hitEnergy;
+        }
+        sigma2ToLine += TVector3(toAddX, toAddY, toAddZ);
+    }
+
+    if (ponderateByEnergy) {
+        auto trackEnergy = track->GetEnergy();
+        sigma2ToLine *= 1.0 / trackEnergy;
+    } else {
+        sigma2ToLine *= 1.0 / track->GetVolumeHits()->GetNumberOfHits();
+    }
+    TVector3 sigmaToLine =
+        TVector3(TMath::Sqrt(sigma2ToLine.X()), TMath::Sqrt(sigma2ToLine.Y()), TMath::Sqrt(sigma2ToLine.Z()));
+    return sigmaToLine;
+}