Merge pull request #31 from rest-for-physics/trackAna

Adding new observables to TRestTrackAnalysisProcess and bug fix on TRestTrackLineAnalysisProcess
rest-for-physics · Mar 28, 2023 · d1c43ce · d1c43ce
2 parents 9c63e85 + e820bb0
commit d1c43ce
Show file tree

Hide file tree

Showing 2 changed files with 18 additions and 3 deletions.
diff --git a/src/TRestTrackAnalysisProcess.cxx b/src/TRestTrackAnalysisProcess.cxx
@@ -872,10 +872,13 @@ TRestEvent* TRestTrackAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
     Double_t tckMaxXZ_SigmaZ = 0;
     Double_t tckMaxYZ_SigmaY = 0;
     Double_t tckMaxYZ_SigmaZ = 0;
+    Double_t tckMaxXZ_nHits = 0;
+    Double_t tckMaxYZ_nHits = 0;
     Double_t tckMaxXYZ_gausSigmaX = 0, tckMaxXYZ_gausSigmaY = 0, tckMaxXYZ_gausSigmaZ = 0;
     Double_t tckMaxXZ_gausSigmaX = 0, tckMaxXZ_gausSigmaZ_XZ = 0;
     Double_t tckMaxYZ_gausSigmaY = 0, tckMaxYZ_gausSigmaZ_YZ = 0;
     Double_t tckMaxTrack_XYZ_GaussSigmaZ = 0;
+    Double_t tckMaxTrack_XYZ_SigmaZ2 = 0;
     Double_t tckMaxTrack_XYZ_skewXY = 0, tckMaxTrack_XYZ_skewZ = 0;
 
     if (fInputTrackEvent->GetMaxEnergyTrack()) {
@@ -903,6 +906,7 @@ TRestEvent* TRestTrackAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
         tckMaxXZ_SigmaZ = fInputTrackEvent->GetMaxEnergyTrack("X")->GetHits()->GetSigmaZ2();
         tckMaxXZ_gausSigmaX = fInputTrackEvent->GetMaxEnergyTrack("X")->GetHits()->GetGaussSigmaX();
         tckMaxXZ_gausSigmaZ_XZ = fInputTrackEvent->GetMaxEnergyTrack("X")->GetHits()->GetGaussSigmaZ();
+        tckMaxXZ_nHits = fInputTrackEvent->GetMaxEnergyTrack("X")->GetNumberOfHits();
         RESTDebug << "id: " << fInputTrackEvent->GetID() << " " << fInputTrackEvent->GetSubEventTag()
                   << " tckMaxEnX: " << tckMaxEnX << RESTendl;
     }
@@ -912,13 +916,15 @@ TRestEvent* TRestTrackAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
     SetObservableValue((string) "MaxTrack_XZ_SigmaZ", tckMaxXZ_SigmaZ);
     SetObservableValue((string) "MaxTrack_XZ_GaussSigmaX", tckMaxXZ_gausSigmaX);
     SetObservableValue((string) "MaxTrack_XZ_GaussSigmaZ", tckMaxXZ_gausSigmaZ_XZ);
+    SetObservableValue((string) "MaxTrack_XZ_nHits", tckMaxXZ_nHits);
 
     if (fInputTrackEvent->GetMaxEnergyTrack("Y")) {
         tckMaxEnY = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetEnergy();
         tckMaxYZ_SigmaY = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetHits()->GetSigmaY();
         tckMaxYZ_SigmaZ = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetHits()->GetSigmaZ2();
         tckMaxYZ_gausSigmaY = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetHits()->GetGaussSigmaY();
         tckMaxYZ_gausSigmaZ_YZ = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetHits()->GetGaussSigmaZ();
+        tckMaxYZ_nHits = fInputTrackEvent->GetMaxEnergyTrack("Y")->GetNumberOfHits();
         RESTDebug << "id: " << fInputTrackEvent->GetID() << " " << fInputTrackEvent->GetSubEventTag()
                   << " tckMaxEnY: " << tckMaxEnY << RESTendl;
     }
@@ -928,6 +934,7 @@ TRestEvent* TRestTrackAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
     SetObservableValue((string) "MaxTrack_YZ_SigmaZ", tckMaxYZ_SigmaZ);
     SetObservableValue((string) "MaxTrack_YZ_GaussSigmaY", tckMaxYZ_gausSigmaY);
     SetObservableValue((string) "MaxTrack_YZ_GaussSigmaZ", tckMaxYZ_gausSigmaZ_YZ);
+    SetObservableValue((string) "MaxTrack_YZ_nHits", tckMaxYZ_nHits);
 
     SetObservableValue("MaxTrackxySigmaGausBalance", (tckMaxXZ_gausSigmaX - tckMaxYZ_gausSigmaY) /
                                                          (tckMaxXZ_gausSigmaX + tckMaxYZ_gausSigmaY));
@@ -981,9 +988,11 @@ TRestEvent* TRestTrackAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
         }
     }
     tckMaxTrack_XYZ_GaussSigmaZ = hits.GetGaussSigmaZ();
+    tckMaxTrack_XYZ_SigmaZ2 = hits.GetSigmaZ2();
     tckMaxTrack_XYZ_skewZ = hits.GetSkewZ();
     tckMaxTrack_XYZ_skewXY = hits.GetSkewXY();
     SetObservableValue((string) "MaxTrack_XYZ_GaussSigmaZ", tckMaxTrack_XYZ_GaussSigmaZ);
+    SetObservableValue((string) "MaxTrack_XYZ_SigmaZ2", tckMaxTrack_XYZ_SigmaZ2);
     SetObservableValue((string) "MaxTrack_XYZ_skewZ", tckMaxTrack_XYZ_skewZ);
     SetObservableValue((string) "MaxTrack_XYZ_skewXY", tckMaxTrack_XYZ_skewXY);
     /* }}} */

diff --git a/src/TRestTrackLineAnalysisProcess.cxx b/src/TRestTrackLineAnalysisProcess.cxx
@@ -39,6 +39,7 @@
 /// ### Observables
 /// * **trackBalanceXZ**: Track balance between the most energetic track and all tracks in the XZ projection
 /// * **trackBalanceYZ**: Track balance between the most energetic track and all tracks in the YZ projection
+/// * **trackBalance**: Total track balance between the most energetic tracks (XZ + YZ) and all tracks
 /// * **originX**: Track origin in the X cordinate
 /// * **originY**: Track origin in the Y cordinate
 /// * **originZ**: Track origin in the Z cordinate
@@ -56,6 +57,7 @@
 ///            verboseLevel="silent" >
 ///                   <observable name="trackBalanceXZ" value="ON" />
 ///                   <observable name="trackBalanceYZ" value="ON" />
+///                   <observable name="trackBalance" value="ON" />
 ///                   <observable name="originX" value="ON" />
 ///                   <observable name="originY" value="ON" />
 ///                   <observable name="originZ" value="ON" />
@@ -138,7 +140,7 @@ TRestEvent* TRestTrackLineAnalysisProcess::ProcessEvent(TRestEvent* inputEvent)
     Double_t angle = -10;
     bool downwards = true;
     Double_t trackEnergyX = 0, trackEnergyY = 0;
-    Double_t trackBalanceX = 0, trackBalanceY = 0;
+    Double_t trackBalanceX = 0, trackBalanceY = 0, trackBalance = 0;
 
     if (tckX && tckY) {
         // Retreive origin and end of the track for the XZ projection
@@ -160,15 +162,19 @@ TRestEvent* TRestTrackLineAnalysisProcess::ProcessEvent(TRestEvent* inputEvent)
         trackEnergyX = tckX->GetEnergy();
         trackEnergyY = tckY->GetEnergy();
 
-        if (trackEnergyX > 0) trackBalanceX = fTrackEvent->GetEnergy("X") / trackEnergyX;
-        if (trackEnergyY > 0) trackBalanceY = fTrackEvent->GetEnergy("Y") / trackEnergyY;
+        if (trackEnergyX > 0) trackBalanceX = trackEnergyX / fTrackEvent->GetEnergy("X");
+        if (trackEnergyY > 0) trackBalanceY = trackEnergyY / fTrackEvent->GetEnergy("Y");
+        if (trackEnergyX > 0 && trackEnergyY > 0)
+            trackBalance =
+                (trackEnergyX + trackEnergyY) / (fTrackEvent->GetEnergy("X") + fTrackEvent->GetEnergy("Y"));
     }
 
     Double_t trackEnergy = trackEnergyX + trackEnergyY;
 
     // A new value for each observable is added
     SetObservableValue("trackBalanceX", trackBalanceX);
     SetObservableValue("trackBalanceY", trackBalanceY);
+    SetObservableValue("trackBalance", trackBalance);
     SetObservableValue("originX", orig.X());
     SetObservableValue("originY", orig.Y());
     SetObservableValue("originZ", orig.Z());