Merge pull request #104 from rest-for-physics/lobis-readout-remove-hits

TRestDetectorHitsReadoutAnalysisProcess: Work correctly with multiple hit types
rest-for-physics · May 7, 2024 · 3ca09f8 · 3ca09f8
2 parents 985cd11 + ab93f20
commit 3ca09f8
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 6 deletions.
diff --git a/inc/TRestDetectorHitsReadoutAnalysisProcess.h b/inc/TRestDetectorHitsReadoutAnalysisProcess.h
@@ -16,8 +16,8 @@
 //! An analysis REST process to extract valuable information from Hits type of data.
 class TRestDetectorHitsReadoutAnalysisProcess : public TRestEventProcess {
    private:
-    TRestDetectorHitsEvent* fInputHitsEvent;   //!
-    TRestDetectorHitsEvent* fOutputHitsEvent;  //!
+    TRestDetectorHitsEvent* fInputHitsEvent = nullptr;   //!
+    TRestDetectorHitsEvent* fOutputHitsEvent = nullptr;  //!
 
     void InitFromConfigFile() override;
     void Initialize() override;

diff --git a/src/TRestDetectorHitsReadoutAnalysisProcess.cxx b/src/TRestDetectorHitsReadoutAnalysisProcess.cxx
@@ -15,12 +15,11 @@ TRestEvent* TRestDetectorHitsReadoutAnalysisProcess::ProcessEvent(TRestEvent* in
     vector<double> hitEnergy;
     double energyInFiducial = 0;
 
-    for (size_t hitIndex = 0; hitIndex < fInputHitsEvent->GetNumberOfHits(); hitIndex++) {
+    for (int hitIndex = 0; hitIndex < static_cast<int>(fInputHitsEvent->GetNumberOfHits()); hitIndex++) {
         const auto position = fInputHitsEvent->GetPosition(hitIndex);
         const auto energy = fInputHitsEvent->GetEnergy(hitIndex);
         const auto time = fInputHitsEvent->GetTime(hitIndex);
         const auto type = fInputHitsEvent->GetType(hitIndex);
-        fOutputHitsEvent->AddHit(position, energy, time, type);
 
         if (energy <= 0) {
             // this should never happen
@@ -31,8 +30,12 @@ TRestEvent* TRestDetectorHitsReadoutAnalysisProcess::ProcessEvent(TRestEvent* in
 
         const auto daqId = fReadout->GetDaqId(position, false);
         const auto channelType = fReadout->GetTypeForChannelDaqId(daqId);
+        const bool isValidHit = channelType == fChannelType;
+
+        // we need to add all hits to preserve the input event
+        fOutputHitsEvent->AddHit(position, energy, time, type);
 
-        if (channelType != fChannelType) {
+        if (!isValidHit) {
             continue;
         }
 
@@ -53,6 +56,7 @@ TRestEvent* TRestDetectorHitsReadoutAnalysisProcess::ProcessEvent(TRestEvent* in
     const double readoutEnergy = accumulate(hitEnergy.begin(), hitEnergy.end(), 0.0);
 
     if (fRemoveZeroEnergyEvents && readoutEnergy <= 0) {
+        // events not depositing any energy in the readout are removed
         return nullptr;
     }
 
@@ -76,6 +80,26 @@ TRestEvent* TRestDetectorHitsReadoutAnalysisProcess::ProcessEvent(TRestEvent* in
     positionSigma.SetY(sqrt(positionSigma.Y()));
     positionSigma.SetZ(sqrt(positionSigma.Z()));
 
+    const auto positionSigmaXY =
+        sqrt(positionSigma.X() * positionSigma.X() + positionSigma.Y() * positionSigma.Y());
+    const auto positionSigmaXYBalance =
+        (positionSigma.X() - positionSigma.Y()) / (positionSigma.X() + positionSigma.Y());
+
+    TVector3 positionSkewness;
+    for (size_t i = 0; i < hitPosition.size(); i++) {
+        TVector3 diff3 = hitPosition[i] - positionAverage;
+        diff3.SetX(diff3.X() * diff3.X() * diff3.X());
+        diff3.SetY(diff3.Y() * diff3.Y() * diff3.Y());
+        diff3.SetZ(diff3.Z() * diff3.Z() * diff3.Z());
+        positionSkewness += diff3 * hitEnergy[i];
+    }
+    positionSkewness *= 1.0 / readoutEnergy;
+    const auto positionSkewnessXY =
+        (positionSkewness.X() + positionSkewness.Y()) / (positionSigmaXY * positionSigmaXY * positionSigmaXY);
+    positionSkewness.SetX(positionSkewness.X() / (positionSigma.X() * positionSigma.X() * positionSigma.X()));
+    positionSkewness.SetY(positionSkewness.Y() / (positionSigma.Y() * positionSigma.Y() * positionSigma.Y()));
+    positionSkewness.SetZ(positionSkewness.Z() / (positionSigma.Z() * positionSigma.Z() * positionSigma.Z()));
+
     SetObservableValue("readoutEnergy", readoutEnergy);
     SetObservableValue("readoutNumberOfHits", hitEnergy.size());
 
@@ -86,6 +110,14 @@ TRestEvent* TRestDetectorHitsReadoutAnalysisProcess::ProcessEvent(TRestEvent* in
     SetObservableValue("readoutSigmaPositionX", positionSigma.X());
     SetObservableValue("readoutSigmaPositionY", positionSigma.Y());
     SetObservableValue("readoutSigmaPositionZ", positionSigma.Z());
+    SetObservableValue("readoutSigmaPositionXY", positionSigmaXY);
+
+    SetObservableValue("readoutSigmaPositionXYBalance", positionSigmaXYBalance);
+
+    SetObservableValue("readoutSkewnessPositionX", positionSkewness.X());
+    SetObservableValue("readoutSkewnessPositionY", positionSkewness.Y());
+    SetObservableValue("readoutSkewnessPositionZ", positionSkewness.Z());
+    SetObservableValue("readoutSkewnessPositionXY", positionSkewnessXY);
 
     SetObservableValue("readoutEnergyInFiducial", energyInFiducial);
 
@@ -100,7 +132,7 @@ void TRestDetectorHitsReadoutAnalysisProcess::InitProcess() {
         exit(1);
     }
 
-    if (fChannelType == "") {
+    if (fChannelType.empty()) {
         cerr << "TRestDetectorHitsReadoutAnalysisProcess::InitProcess() : "
              << "Channel type not defined" << endl;
         exit(1);