TRestDetectorHitmapAnalysisProcess.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 "TRestDetectorHitmapAnalysisProcess.h"
 
using namespace std;
 
ClassImp(TRestDetectorHitmapAnalysisProcess);
 
TRestDetectorHitmapAnalysisProcess::TRestDetectorHitmapAnalysisProcess() { Initialize(); }
 
TRestDetectorHitmapAnalysisProcess::~TRestDetectorHitmapAnalysisProcess() {}
 
TRestDetectorHitmapAnalysisProcess::TRestDetectorHitmapAnalysisProcess(const char* configFilename) {
    Initialize();
    LoadConfigFromFile(configFilename);
}
 
void TRestDetectorHitmapAnalysisProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fHitsEvent = nullptr;
}
 
TRestEvent* TRestDetectorHitmapAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
    fHitsEvent = (TRestDetectorHitsEvent*)inputEvent;
 
    TVector3 hitMean = fHitsEvent->GetMeanPosition();
 
    TVector3 transformedMean = hitMean;
 
    for (const auto& t : fTransformations) {
        HitTransformation tr;
        tr = GetTransformation(t);
 
        if (tr.type == "specular") transformedMean = Specular(transformedMean, tr);
        if (tr.type == "rotation") transformedMean = Rotation(transformedMean, tr);
        if (tr.type == "translation") transformedMean = Translation(transformedMean, tr);
    }
 
    SetObservableValue("xMean", transformedMean.X());
    SetObservableValue("yMean", transformedMean.Y());
    SetObservableValue("zMean", transformedMean.Z());
 
    return fHitsEvent;
}
 
TVector3 TRestDetectorHitmapAnalysisProcess::Specular(const TVector3& pos, const HitTransformation& tr) {
    if (tr.type != "specular") return {0, 0, 0};
 
    TVector3 V = pos - tr.position;
    Double_t vByN = V.Dot(tr.vector);
    TVector3 reflectionVector = V - 2 * vByN * tr.vector;
    return tr.position + reflectionVector;
}
 
TVector3 TRestDetectorHitmapAnalysisProcess::Rotation(const TVector3& pos, const HitTransformation& tr) {
    if (tr.type != "rotation") return {0, 0, 0};
 
    TVector3 position = pos - tr.position;
    position.Rotate(tr.angle, tr.vector);
    position += tr.position;
 
    return position;
}
 
TVector3 TRestDetectorHitmapAnalysisProcess::Translation(const TVector3& pos, const HitTransformation& tr) {
    if (tr.type != "translation") return {0, 0, 0};
 
    TVector3 position = pos + tr.vector;
    return position;
}
 
HitTransformation TRestDetectorHitmapAnalysisProcess::GetTransformation(const std::string& name) {
    for (const auto& t : fTransDefinitions) {
        if (t.name == name) return t;
    }
 
    HitTransformation Tr;
    return Tr;
}
 
void TRestDetectorHitmapAnalysisProcess::InitFromConfigFile() {
    TRestEventProcess::InitFromConfigFile();
 
    auto specularDef = GetElement("specular");
    while (specularDef) {
        HitTransformation trans;
        trans.type = "specular";
 
        string name = GetFieldValue("name", specularDef);
        if (name == "Not defined") {
            RESTError << "The `name` field was not found in specular definition!" << RESTendl;
            exit(1);
        }
        trans.name = name;
 
        TVector3 position = Get3DVectorParameterWithUnits("position", specularDef);
        if (position == TVector3(-1, -1, -1)) {
            RESTError << "The `position` field was not found in specular definition!" << RESTendl;
            exit(1);
        }
        trans.position = position;
 
        TVector3 vect = Get3DVectorParameterWithUnits("vector", specularDef);
        if (vect == TVector3(-1, -1, -1)) {
            RESTError << "The `vector` field was not found in specular definition!" << RESTendl;
            exit(1);
        }
        trans.vector = vect;
 
        fTransDefinitions.push_back(trans);
        specularDef = GetNextElement(specularDef);
    }
 
    auto rotationDef = GetElement("rotation");
    while (rotationDef) {
        HitTransformation trans;
        trans.type = "rotation";
 
        string name = GetFieldValue("name", rotationDef);
        if (name == "Not defined") {
            RESTError << "The `name` field was not found in rotation definition!" << RESTendl;
            exit(1);
        }
        trans.name = name;
 
        TVector3 position = Get3DVectorParameterWithUnits("position", rotationDef);
        if (position == TVector3(-1, -1, -1)) {
            RESTError << "The `position` field was not found in rotation definition!" << RESTendl;
            exit(1);
        }
        trans.position = position;
 
        TVector3 vect = Get3DVectorParameterWithUnits("vector", rotationDef);
        if (vect == TVector3(-1, -1, -1)) {
            RESTError << "The `vector` field was not found in rotation definition!" << RESTendl;
            exit(1);
        }
        trans.vector = vect;
 
        Double_t angle = GetDblParameterWithUnits("angle", rotationDef);
        if (angle == PARAMETER_NOT_FOUND_DBL) {
            RESTError << "The `angle` field was not found in rotation definition!" << RESTendl;
            exit(1);
        }
        trans.angle = angle;
 
        fTransDefinitions.push_back(trans);
        rotationDef = GetNextElement(rotationDef);
    }
 
    auto translationDef = GetElement("translation");
    while (translationDef) {
        HitTransformation trans;
        trans.type = "translation";
 
        string name = GetFieldValue("name", translationDef);
        if (name == "Not defined") {
            RESTError << "The `name` field was not found in translation definition!" << RESTendl;
            exit(1);
        }
        trans.name = name;
 
        TVector3 vect = Get3DVectorParameterWithUnits("vector", translationDef);
        if (vect == TVector3(-1, -1, -1)) {
            RESTError << "The `vector` field was not found in translation definition!" << RESTendl;
            exit(1);
        }
        trans.vector = vect;
 
        fTransDefinitions.push_back(trans);
        translationDef = GetNextElement(translationDef);
    }
}
 
void TRestDetectorHitmapAnalysisProcess::PrintMetadata() {
    BeginPrintProcess();
 
    if (!fTransDefinitions.empty()) {
        RESTMetadata << " List of transformation definitions: " << RESTendl;
        RESTMetadata << " ----------------------------------- " << RESTendl;
        RESTMetadata << " " << RESTendl;
    } else
        RESTMetadata << " No transformations found! " << RESTendl;
 
    for (const auto& t : fTransDefinitions) {
        if (t.type == "specular") {
            RESTMetadata << " - Specular transformation: " << RESTendl;
            RESTMetadata << "   + Plane position : (" << t.position.X() << ", " << t.position.Y() << ", "
                         << t.position.Z() << ") mm" << RESTendl;
            RESTMetadata << "   + Plane normal : (" << t.vector.X() << ", " << t.vector.Y() << ", "
                         << t.vector.Z() << ")" << RESTendl;
            RESTMetadata << " " << RESTendl;
        }
 
        if (t.type == "rotation") {
            RESTMetadata << " - Rotation transformation: " << RESTendl;
            RESTMetadata << "   + Rotation center : (" << t.position.X() << ", " << t.position.Y() << ", "
                         << t.position.Z() << ") mm" << RESTendl;
            RESTMetadata << "   + Rotation axis : (" << t.vector.X() << ", " << t.vector.Y() << ", "
                         << t.vector.Z() << ")" << RESTendl;
            RESTMetadata << "   + Rotation angle : " << t.angle * units("degrees") << " degrees" << RESTendl;
            RESTMetadata << " " << RESTendl;
        }
 
        if (t.type == "translation") {
            RESTMetadata << " - Translation transformation: " << RESTendl;
            RESTMetadata << "   + Translation vector : (" << t.vector.X() << ", " << t.vector.Y() << ", "
                         << t.vector.Z() << ") mm" << RESTendl;
            RESTMetadata << " " << RESTendl;
        }
    }
 
    RESTMetadata << " Order of transformations: " << RESTendl;
    RESTMetadata << " ------------------------- " << RESTendl;
    RESTMetadata << " " << RESTendl;
 
    int n = 0;
    for (const auto& t : fTransformations) {
        n++;
        RESTMetadata << n << ". " << t << RESTendl;
    }
 
    EndPrintProcess();
}