TRestTrack3DAnalysisProcess.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 "TRestTrack3DAnalysisProcess.h"
using namespace std;
 
// Comparator function for sorting in descending order
bool sortByValueDescending3D(const std::pair<int, Double_t>& a, const std::pair<int, Double_t>& b) {
    return a.second > b.second;  // Change to < for ascending order
}
 
ClassImp(TRestTrack3DAnalysisProcess);
 
TRestTrack3DAnalysisProcess::TRestTrack3DAnalysisProcess() { Initialize(); }
 
TRestTrack3DAnalysisProcess::TRestTrack3DAnalysisProcess(const char* configFilename) {
    Initialize();
    if (LoadConfigFromFile(configFilename)) LoadDefaultConfig();
}
 
TRestTrack3DAnalysisProcess::~TRestTrack3DAnalysisProcess() { delete fTrackEvent; }
 
void TRestTrack3DAnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); }
 
void TRestTrack3DAnalysisProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fTrackEvent = new TRestTrackEvent();
}
 
void TRestTrack3DAnalysisProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
}
 
void TRestTrack3DAnalysisProcess::InitProcess() {}
 
TRestEvent* TRestTrack3DAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
    fTrackEvent = (TRestTrackEvent*)inputEvent;
 
 
    map<int, Double_t> XYZ_Energy;
 
    map<int, int> XYZ_NHits;
 
    map<int, Double_t> XYZ_SigmaX;
    map<int, Double_t> XYZ_SigmaY;
    map<int, Double_t> XYZ_SigmaZ;
 
    map<int, Double_t> XYZ_SkewXY;
    map<int, Double_t> XYZ_SkewZ;
 
    map<int, Double_t> XYZ_GaussSigmaX;
    map<int, Double_t> XYZ_GaussSigmaY;
    map<int, Double_t> XYZ_GaussSigmaZ;
 
    map<int, Double_t> XYZ_Length;
 
    map<int, Double_t> XYZ_Volume;
 
    map<int, Double_t> XYZ_MeanX;
    map<int, Double_t> XYZ_MeanY;
    map<int, Double_t> XYZ_MeanZ;
 
    Double_t XYZ_FirstSecondTracksDistance;
 
 
    for (int tck = 0; tck < fTrackEvent->GetNumberOfTracks(); tck++) {
        if (!fTrackEvent->isTopLevel(tck)) continue;
 
        TRestTrack* t = fTrackEvent->GetTrack(tck);
 
        if (t->isXYZ()) {
            XYZ_NHits[t->GetTrackID()] = t->GetNumberOfHits();
            XYZ_Energy[t->GetTrackID()] = t->GetTrackEnergy();
            XYZ_SigmaX[t->GetTrackID()] = t->GetHits()->GetSigmaX();
            XYZ_SigmaY[t->GetTrackID()] = t->GetHits()->GetSigmaY();
            XYZ_SigmaZ[t->GetTrackID()] = t->GetHits()->GetSigmaZ2();
            XYZ_GaussSigmaX[t->GetTrackID()] = t->GetHits()->GetGaussSigmaX();
            XYZ_GaussSigmaY[t->GetTrackID()] = t->GetHits()->GetGaussSigmaY();
            XYZ_GaussSigmaZ[t->GetTrackID()] = t->GetHits()->GetGaussSigmaZ();
            XYZ_Length[t->GetTrackID()] = t->GetLength();
            XYZ_Volume[t->GetTrackID()] = t->GetVolume();
            XYZ_MeanX[t->GetTrackID()] = t->GetMeanPosition().X();
            XYZ_MeanY[t->GetTrackID()] = t->GetMeanPosition().Y();
            XYZ_MeanZ[t->GetTrackID()] = t->GetMeanPosition().Z();
            XYZ_SkewXY[t->GetTrackID()] = t->GetHits()->GetSkewXY();
            XYZ_SkewZ[t->GetTrackID()] = t->GetHits()->GetSkewZ();
        } else {
            XYZ_NHits[t->GetTrackID()] = 0;
            XYZ_Energy[t->GetTrackID()] = 0;
            XYZ_SigmaX[t->GetTrackID()] = 0;
            XYZ_SigmaY[t->GetTrackID()] = 0;
            XYZ_SigmaZ[t->GetTrackID()] = 0;
            XYZ_GaussSigmaX[t->GetTrackID()] = 0;
            XYZ_GaussSigmaY[t->GetTrackID()] = 0;
            XYZ_GaussSigmaZ[t->GetTrackID()] = 0;
            XYZ_Length[t->GetTrackID()] = 0;
            XYZ_Volume[t->GetTrackID()] = 0;
            XYZ_MeanX[t->GetTrackID()] = 0;
            XYZ_MeanY[t->GetTrackID()] = 0;
            XYZ_MeanZ[t->GetTrackID()] = 0;
            XYZ_SkewXY[t->GetTrackID()] = 0;
            XYZ_SkewZ[t->GetTrackID()] = 0;
        }
    }
 
    vector<pair<int, Double_t>> energies;
 
    map<int, Double_t>::iterator it;
    for (it = XYZ_Energy.begin(); it != XYZ_Energy.end(); it++) {
        energies.emplace_back(make_pair(it->first, it->second));
    }
 
    sort(energies.begin(), energies.end(), sortByValueDescending3D);
 
    if (fTrackEvent->GetNumberOfTracks() > 1) {
        Double_t dX = 0, dY = 0, dZ = 0;
 
        dX = abs(XYZ_MeanX[energies[0].first] - XYZ_MeanX[energies[1].first]);
        dY = abs(XYZ_MeanY[energies[0].first] - XYZ_MeanY[energies[1].first]);
        dZ = abs(XYZ_MeanZ[energies[0].first] - XYZ_MeanZ[energies[1].first]);
 
        XYZ_FirstSecondTracksDistance = TMath::Sqrt(dX * dX + dY * dY + dZ * dZ);
    } else {
        XYZ_FirstSecondTracksDistance = 0;
    }
 
 
    // --- Number of tracks and energy --- //
    SetObservableValue("NTracks", fTrackEvent->GetNumberOfTracks());
    SetObservableValue("TotalEnergy", fTrackEvent->GetEnergy());
 
    // --- Map observables --- //
    SetObservableValue("Map_XYZ_NHits", XYZ_NHits);
    SetObservableValue("Map_XYZ_Energy", XYZ_Energy);
    SetObservableValue("Map_XYZ_SigmaX", XYZ_SigmaX);
    SetObservableValue("Map_XYZ_SigmaY", XYZ_SigmaY);
    SetObservableValue("Map_XYZ_SigmaZ", XYZ_SigmaZ);
    SetObservableValue("Map_XYZ_GaussSigmaX", XYZ_GaussSigmaX);
    SetObservableValue("Map_XYZ_GaussSigmaY", XYZ_GaussSigmaY);
    SetObservableValue("Map_XYZ_GaussSigmaZ", XYZ_GaussSigmaZ);
    SetObservableValue("Map_XYZ_Length", XYZ_Length);
    SetObservableValue("Map_XYZ_Volume", XYZ_Volume);
    SetObservableValue("Map_XYZ_MeanX", XYZ_MeanX);
    SetObservableValue("Map_XYZ_MeanY", XYZ_MeanY);
    SetObservableValue("Map_XYZ_MeanZ", XYZ_MeanZ);
    SetObservableValue("Map_XYZ_SkewXY", XYZ_SkewXY);
    SetObservableValue("Map_XYZ_SkewZ", XYZ_SkewZ);
 
    // Copy the MaxTrack keys immediately after checking the vector
    if (!energies.empty()) {
        // --- Max track observables --- //
        int energies0FirstKey =
            energies[0].first;  // Declare Keys outside to avoid error when accessing "energies[0].first"...
        Double_t energies0SecondKey = energies[0].second;
 
        SetObservableValue("MaxTrack_XYZ_NHits", XYZ_NHits[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_Energy", XYZ_Energy[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_SigmaX", XYZ_SigmaX[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_SigmaY", XYZ_SigmaY[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_SigmaZ", XYZ_SigmaZ[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_GaussSigmaX", XYZ_GaussSigmaX[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_GaussSigmaY", XYZ_GaussSigmaY[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_GaussSigmaZ", XYZ_GaussSigmaZ[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_Length", XYZ_Length[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_Volume", XYZ_Volume[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_MeanX", XYZ_MeanX[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_MeanY", XYZ_MeanY[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_MeanZ", XYZ_MeanZ[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_SkewZ", XYZ_SkewXY[energies0FirstKey]);
        SetObservableValue("MaxTrack_XYZ_SkewZ", XYZ_SkewZ[energies0FirstKey]);
 
        SetObservableValue("MaxTrack_XYZ_MaxTrackEnergyPercentage",
                           (energies0SecondKey) / fTrackEvent->GetEnergy());
 
        // --- Second max track observables --- //
        int energies1FirstKey = energies[1].first;
 
        SetObservableValue("SecondMaxTrack_XYZ_NHits", XYZ_NHits[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_Energy", XYZ_Energy[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_SigmaX", XYZ_SigmaX[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_SigmaY", XYZ_SigmaY[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_SigmaZ", XYZ_SigmaZ[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_GaussSigmaX", XYZ_GaussSigmaX[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_GaussSigmaY", XYZ_GaussSigmaY[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_GaussSigmaZ", XYZ_GaussSigmaZ[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_Length", XYZ_Length[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_Volume", XYZ_Volume[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_MeanX", XYZ_MeanX[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_MeanY", XYZ_MeanY[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_MeanZ", XYZ_MeanZ[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_SkewZ", XYZ_SkewXY[energies1FirstKey]);
        SetObservableValue("SecondMaxTrack_XYZ_SkewZ", XYZ_SkewZ[energies1FirstKey]);
 
        Double_t energies1SecondKey = energies[1].second;
 
        if (fTrackEvent->GetNumberOfTracks() > 2) {
            SetObservableValue("SecondMaxTrack_XYZ_MaxTrackEnergyPercentage",
                               (energies1SecondKey) / fTrackEvent->GetEnergy());
        } else {
            SetObservableValue("SecondMaxTrack_XYZ_MaxTrackEnergyPercentage", 0.0);
        }
    } else {
        std::cerr << "Error: energies is empty. Some observables will not be set." << std::endl;
    }
 
    // --- Distance observables between first two tracks --- //
    SetObservableValue("XYZ_FirstSecondTracksDistance", XYZ_FirstSecondTracksDistance);
 
    return fTrackEvent;
}
 
void TRestTrack3DAnalysisProcess::EndProcess() {}
 
void TRestTrack3DAnalysisProcess::InitFromConfigFile() {}