TRestDetectorHitsAnalysisProcess.cxx

/*************************************************************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see http://gifna.unizar.es/trex                  *
 *                                                                       *
 * REST is free software: you can redistribute it and/or modify          *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * REST is distributed in the hope that it will be useful,               *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have a copy of the GNU General Public License along with   *
 * REST in $REST_PATH/LICENSE.                                           *
 * If not, see http://www.gnu.org/licenses/.                             *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
 
#include "TRestDetectorHitsAnalysisProcess.h"
 
using namespace std;
 
ClassImp(TRestDetectorHitsAnalysisProcess);
 
TRestDetectorHitsAnalysisProcess::TRestDetectorHitsAnalysisProcess() { Initialize(); }
 
TRestDetectorHitsAnalysisProcess::TRestDetectorHitsAnalysisProcess(const char* configFilename) {
    Initialize();
    if (LoadConfigFromFile(configFilename) == -1) {
        LoadDefaultConfig();
    }
}
 
TRestDetectorHitsAnalysisProcess::~TRestDetectorHitsAnalysisProcess() { delete fOutputHitsEvent; }
 
void TRestDetectorHitsAnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); }
 
void TRestDetectorHitsAnalysisProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fInputHitsEvent = nullptr;
    fOutputHitsEvent = new TRestDetectorHitsEvent();
 
    fPrismFiducial = false;
    fCylinderFiducial = false;
}
 
void TRestDetectorHitsAnalysisProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name) == -1) LoadDefaultConfig();
}
 
void TRestDetectorHitsAnalysisProcess::InitProcess() { TRestEventProcess::ReadObservables(); }
 
TRestEvent* TRestDetectorHitsAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
    fInputHitsEvent = (TRestDetectorHitsEvent*)inputEvent;
 
    TString obsName;
 
    TRestHits* hits = fInputHitsEvent->GetHits();
    for (unsigned int n = 0; n < hits->GetNumberOfHits(); n++) {
        Double_t eDep = hits->GetEnergy(n);
 
        Double_t x = hits->GetX(n);
        Double_t y = hits->GetY(n);
        Double_t z = hits->GetZ(n);
 
        auto time = hits->GetTime(n);
        auto type = hits->GetType(n);
 
        fOutputHitsEvent->AddHit(x, y, z, eDep, time, type);
    }
 
    Double_t energy = fOutputHitsEvent->GetEnergy();
    TVector3 meanPosition = fOutputHitsEvent->GetMeanPosition();
    Double_t sigmaX = fOutputHitsEvent->GetSigmaX();
    Double_t sigmaY = fOutputHitsEvent->GetSigmaY();
    Double_t sigmaXY2 = fOutputHitsEvent->GetSigmaXY2();
    Double_t sigmaZ2 = fOutputHitsEvent->GetSigmaZ2();
    Double_t skewXY = fOutputHitsEvent->GetSkewXY();
    Double_t skewZ = fOutputHitsEvent->GetSkewZ();
    Double_t energyX = fOutputHitsEvent->GetEnergyX();
    Double_t energyY = fOutputHitsEvent->GetEnergyY();
    Double_t maxEnergy = fOutputHitsEvent->GetMaximumHitEnergy();
    Double_t minEnergy = fOutputHitsEvent->GetMinimumHitEnergy();
    Double_t meanEnergy = fOutputHitsEvent->GetMeanHitEnergy();
    Int_t nHits = fOutputHitsEvent->GetNumberOfHits();
    Int_t nHitsX = fOutputHitsEvent->GetNumberOfHitsX();
    Int_t nHitsY = fOutputHitsEvent->GetNumberOfHitsY();
 
    SetObservableValue("nHits", nHits);
 
    SetObservableValue("nHitsX", nHitsX);
 
    SetObservableValue("nHitsY", nHitsY);
 
    SetObservableValue("balanceXYnHits", (nHitsX - nHitsY) / double(nHitsX + nHitsY));
 
    if ((nHits == nHitsX) || (nHits == nHitsY))
        SetObservableValue("nHitsSizeXY", (double)nHits);
    else
        SetObservableValue("nHitsSizeXY", TMath::Sqrt(nHitsX * nHitsX + nHitsY * nHitsY));
 
    // Checking hits inside fiducial cylinder
    if (fCylinderFiducial) {
        TVector3 meanPositionInCylinder =
            fOutputHitsEvent->GetMeanPositionInCylinder(fFid_x0, fFid_x1, fFid_R);
 
        Int_t isInsideCylinder = 0;
        if (fOutputHitsEvent->anyHitInsideCylinder(fFid_x0, fFid_x1, fFid_R)) isInsideCylinder = 1;
 
        Int_t nCylVol = fOutputHitsEvent->GetNumberOfHitsInsideCylinder(fFid_x0, fFid_x1, fFid_R);
 
        Double_t enCylVol = fOutputHitsEvent->GetEnergyInCylinder(fFid_x0, fFid_x1, fFid_R);
 
        SetObservableValue("isInsideCylindricalVolume", isInsideCylinder);
 
        SetObservableValue("nInsideCylindricalVolume", nCylVol);
 
        SetObservableValue("energyInsideCylindricalVolume", enCylVol);
 
        // mean positions
        SetObservableValue("xMeanInCylinder", meanPositionInCylinder.X());
 
        SetObservableValue("yMeanInCylinder", meanPositionInCylinder.Y());
 
        SetObservableValue("zMeanInCylinder", meanPositionInCylinder.Z());
    }
 
    // Checking hits inside fiducial prism
    if (fPrismFiducial) {
        TVector3 meanPositionInPrism =
            fOutputHitsEvent->GetMeanPositionInPrism(fFid_x0, fFid_x1, fFid_sX, fFid_sY, fFid_theta);
 
        Int_t isInsidePrism = 0;
        if (fOutputHitsEvent->anyHitInsidePrism(fFid_x0, fFid_x1, fFid_sX, fFid_sY, fFid_theta))
            isInsidePrism = 1;
 
        Int_t nPrismVol =
            fOutputHitsEvent->GetNumberOfHitsInsidePrism(fFid_x0, fFid_x1, fFid_sX, fFid_sY, fFid_theta);
 
        Double_t enPrismVol =
            fOutputHitsEvent->GetEnergyInPrism(fFid_x0, fFid_x1, fFid_sX, fFid_sY, fFid_theta);
 
        SetObservableValue("isInsidePrismVolume", isInsidePrism);
 
        SetObservableValue("nInsidePrismVolume", nPrismVol);
 
        SetObservableValue("energyInsidePrismVolume", enPrismVol);
 
        // Mean Positions
 
        SetObservableValue("xMeanInPrism", meanPositionInPrism.X());
 
        SetObservableValue("yMeanInPrism", meanPositionInPrism.Y());
 
        SetObservableValue("zMeanInPrism", meanPositionInPrism.Z());
    }
 
 
    if (fCylinderFiducial) {
        // Adding distances to cylinder wall
        Double_t dToCylWall =
            fOutputHitsEvent->GetClosestHitInsideDistanceToCylinderWall(fFid_x0, fFid_x1, fFid_R);
        Double_t dToCylTop =
            fOutputHitsEvent->GetClosestHitInsideDistanceToCylinderTop(fFid_x0, fFid_x1, fFid_R);
        Double_t dToCylBottom =
            fOutputHitsEvent->GetClosestHitInsideDistanceToCylinderBottom(fFid_x0, fFid_x1, fFid_R);
 
        SetObservableValue("distanceToCylinderWall", dToCylWall);
        SetObservableValue("distanceToCylinderTop", dToCylTop);
        SetObservableValue("distanceToCylinderBottom", dToCylBottom);
    }
 
    if (fPrismFiducial) {
        // Adding distances to prism wall
        Double_t dToPrismWall = fOutputHitsEvent->GetClosestHitInsideDistanceToPrismWall(
            fFid_x0, fFid_x1, fFid_sX, fFid_sY, fFid_theta);
        Double_t dToPrismTop = fOutputHitsEvent->GetClosestHitInsideDistanceToPrismTop(
            fFid_x0, fFid_x1, fFid_sX, fFid_sY, fFid_theta);
        Double_t dToPrismBottom = fOutputHitsEvent->GetClosestHitInsideDistanceToPrismBottom(
            fFid_x0, fFid_x1, fFid_sX, fFid_sY, fFid_theta);
 
        SetObservableValue("distanceToPrismWall", dToPrismWall);
 
        SetObservableValue("distanceToPrismTop", dToPrismTop);
 
        SetObservableValue("distanceToPrismBottom", dToPrismBottom);
    }
 
 
    SetObservableValue("energy", energy);
    SetObservableValue("energyX", energyX);
    SetObservableValue("energyY", energyY);
    SetObservableValue("balanceXYenergy", (energyX - energyY) / (energyX + energyY));
 
    SetObservableValue("maxHitEnergy", maxEnergy);
    SetObservableValue("minHitEnergy", minEnergy);
    SetObservableValue("meanHitEnergy", meanEnergy);
    SetObservableValue("meanHitEnergyBalance", meanEnergy / energy);
 
    SetObservableValue("xMean", meanPosition.X());
 
    SetObservableValue("yMean", meanPosition.Y());
 
    SetObservableValue("zMean", meanPosition.Z());
    SetObservableValue("xy2Sigma", sigmaXY2);
    SetObservableValue("xSigma", sigmaX);
    SetObservableValue("ySigma", sigmaY);
    SetObservableValue("xySigmaBalance", (sigmaX - sigmaY) / (sigmaX + sigmaY));
 
    SetObservableValue("z2Sigma", sigmaZ2);
 
    SetObservableValue("xySkew", skewXY);
    SetObservableValue("zSkew", skewZ);
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
        fOutputHitsEvent->PrintEvent(1000);
        GetChar();
    }
 
    if (fOutputHitsEvent->GetNumberOfHits() == 0) return nullptr;
 
    return fOutputHitsEvent;
}
 
void TRestDetectorHitsAnalysisProcess::EndProcess() {
    // Function to be executed once at the end of the process
    // (after all events have been processed)
 
    // Start by calling the EndProcess function of the abstract class.
    // Comment this if you don't want it.
    // TRestEventProcess::EndProcess();
}
 
void TRestDetectorHitsAnalysisProcess::InitFromConfigFile() {
    fFid_x0 = Get3DVectorParameterWithUnits("fiducial_x0", TVector3(0, 0, 0));
    fFid_x1 = Get3DVectorParameterWithUnits("fiducial_x1", TVector3(0, 0, 0));
 
    fFid_R = GetDblParameterWithUnits("fiducial_R", 1);
    fFid_sX = GetDblParameterWithUnits("fiducial_sX", 1);
    fFid_sY = GetDblParameterWithUnits("fiducial_sY", 1);
    // read angle in degrees
    fFid_theta = StringToDouble(GetParameter("fiducial_theta", "0"));
    // convert it into radians for internal usage
    fFid_theta = std::fmod(fFid_theta, 360) * TMath::DegToRad();
 
    if (GetParameter("cylinderFiducialization", "false") == "true") fCylinderFiducial = true;
 
    if (GetParameter("prismFiducialization", "false") == "true") fPrismFiducial = true;
}