diff --git a/inc/TRestDetectorHitsReadoutAnalysisProcess.h b/inc/TRestDetectorHitsReadoutAnalysisProcess.h index ece71bcd..1aefdc55 100644 --- 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 index 9ba22c46..37fdb49d 100644 --- a/src/TRestDetectorHitsReadoutAnalysisProcess.cxx +++ b/src/TRestDetectorHitsReadoutAnalysisProcess.cxx @@ -15,12 +15,11 @@ TRestEvent* TRestDetectorHitsReadoutAnalysisProcess::ProcessEvent(TRestEvent* in vector hitEnergy; double energyInFiducial = 0; - for (size_t hitIndex = 0; hitIndex < fInputHitsEvent->GetNumberOfHits(); hitIndex++) { + for (int hitIndex = 0; hitIndex < static_cast(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);