src/TRestDetectorReadout.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 "TRestDetectorReadout.h"
 
#include <TFile.h>
 
using namespace std;
 
ClassImp(TRestDetectorReadout);
TRestDetectorReadout::TRestDetectorReadout() { Initialize(); }
 
TRestDetectorReadout::TRestDetectorReadout(const char* configFilename) : TRestMetadata(configFilename) {
    cout << "Loading readout. This might take few seconds" << endl;
    Initialize();
 
    LoadConfigFromFile(fConfigFileName);
}
 
TRestDetectorReadout::TRestDetectorReadout(const char* configFilename, const string& name)
    : TRestMetadata(configFilename) {
    cout << "Loading readout. This might take few seconds" << endl;
    Initialize();
 
    LoadConfigFromFile(fConfigFileName, name);
}
 
void TRestDetectorReadout::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
    fReadoutPlanes.clear();
}
 
TRestDetectorReadout::~TRestDetectorReadout() = default;
 
Int_t TRestDetectorReadout::GetNumberOfModules() {
    Int_t modules = 0;
    for (int p = 0; p < GetNumberOfReadoutPlanes(); p++) {
        modules += fReadoutPlanes[p].GetNumberOfModules();
    }
    return modules;
}
 
Int_t TRestDetectorReadout::GetNumberOfChannels() {
    Int_t channels = 0;
    for (int p = 0; p < GetNumberOfReadoutPlanes(); p++) {
        for (size_t m = 0; m < fReadoutPlanes[p].GetNumberOfModules(); m++) {
            channels += fReadoutPlanes[p][m].GetNumberOfChannels();
        }
    }
    return channels;
}
 
Int_t TRestDetectorReadout::GetModuleDefinitionId(const TString& name) {
    for (unsigned int i = 0; i < fModuleDefinitions.size(); i++) {
        if (fModuleDefinitions[i].GetName() == name) {
            return i;
        }
    }
    return -1;
}
 
TRestDetectorReadoutPlane* TRestDetectorReadout::GetReadoutPlaneWithID(int id) {
    for (int i = 0; i < this->GetNumberOfReadoutPlanes(); i++) {
        if (fReadoutPlanes[i].GetID() == id) {
            return &fReadoutPlanes[i];
        }
    }
 
    return nullptr;
}
 
TRestDetectorReadoutModule* TRestDetectorReadout::GetReadoutModuleWithID(int id) {
    for (int i = 0; i < this->GetNumberOfReadoutPlanes(); i++) {
        TRestDetectorReadoutPlane& plane = fReadoutPlanes[i];
 
        for (size_t j = 0; j < plane.GetNumberOfModules(); j++) {
            if (plane[j].GetModuleID() == id) {
                return &plane[j];
            }
        }
    }
    return nullptr;
}
 
TRestDetectorReadoutChannel* TRestDetectorReadout::GetReadoutChannelWithDaqID(int daqId) {
    for (int p = 0; p < GetNumberOfReadoutPlanes(); p++) {
        for (size_t m = 0; m < fReadoutPlanes[p].GetNumberOfModules(); m++) {
            for (size_t c = 0; c < fReadoutPlanes[p][m].GetNumberOfChannels(); c++) {
                if (fReadoutPlanes[p][m][c].GetDaqID() == daqId) {
                    return &fReadoutPlanes[p][m][c];
                }
            }
        }
    }
    return nullptr;
}
 
TRestDetectorReadoutPlane* TRestDetectorReadout::GetReadoutPlane(int p) {
    if (p < GetNumberOfReadoutPlanes())
        return &fReadoutPlanes[p];
    else {
        RESTWarning << "TRestDetectorReadout::GetReadoutPlane." << RESTendl;
        RESTWarning << "Readout plane index exceeded." << RESTendl;
        RESTWarning << "Index requested : " << p << RESTendl;
        RESTWarning << "Number of readout planes defined : " << GetNumberOfReadoutPlanes() << RESTendl;
    }
 
    return nullptr;
}
 
void TRestDetectorReadout::AddReadoutPlane(const TRestDetectorReadoutPlane& plane) {
    fReadoutPlanes.emplace_back(plane);
}
 
void TRestDetectorReadout::InitFromConfigFile() {
    fMappingNodes = StringToInteger(GetParameter("mappingNodes", "0"));
 
    TiXmlElement* moduleDefinition = GetElement("readoutModule");
    while (moduleDefinition != nullptr) {
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
            cout << "------module-----------------" << endl;
            cout << moduleDefinition << endl;
            cout << "-----------------------------" << endl;
            GetChar();
        }
 
        TRestDetectorReadoutModule module = *ParseModuleDefinition(moduleDefinition);
        module.SetMappingNodes(fMappingNodes);
        module.DoReadoutMapping();
        fModuleDefinitions.push_back(module);
        moduleDefinition = GetNextElement(moduleDefinition);
    }
 
    TiXmlElement* planeDefinition = GetElement("readoutPlane");
    Int_t addedChannels = 0;
    vector<TRestDetectorReadoutModule> moduleVector;
    while (planeDefinition != nullptr) {
        TRestDetectorReadoutPlane plane;
 
        plane.SetID(GetNumberOfReadoutPlanes());
        plane.SetPosition(Get3DVectorParameterWithUnits("position", planeDefinition));
        plane.SetNormal(Get3DVectorParameterWithUnits("normal", planeDefinition));
        plane.SetHeight(GetDblParameterWithUnits("height", planeDefinition));
        plane.SetChargeCollection(StringToDouble(GetFieldValue("chargeCollection", planeDefinition)));
        plane.SetRotation(GetDblParameterWithUnits("rotation", planeDefinition, 0));
        plane.SetName(GetParameter("name", planeDefinition, plane.GetName()));
        plane.SetType(GetParameter("type", planeDefinition, plane.GetType()));
 
        cout << "plane name: " << plane.GetName() << endl;
        cout << "plane type: " << plane.GetType() << endl;
 
        moduleVector.clear();
        TiXmlElement* moduleDefinition = GetElement("addReadoutModule", planeDefinition);
        while (moduleDefinition != nullptr) {
            TString modName = GetFieldValue("name", moduleDefinition);
            Int_t mid = GetModuleDefinitionId(modName);
 
            if (mid == -1) {
                RESTError << "TRestDetectorReadout at <addReadoutModule>. Module name " << modName
                          << " not found!" << RESTendl;
                RESTError << "Please, check spelling" << RESTendl;
                exit(1);
            }
 
            fModuleDefinitions[mid].SetModuleID(StringToInteger(GetFieldValue("id", moduleDefinition)));
            fModuleDefinitions[mid].SetOrigin(StringTo2DVector(GetFieldValue("origin", moduleDefinition)));
            fModuleDefinitions[mid].SetRotation(GetDblParameterWithUnits("rotation", moduleDefinition));
 
            Int_t firstDaqChannel = StringToInteger(GetFieldValue("firstDaqChannel", moduleDefinition));
 
            if (firstDaqChannel == -1) {
                fModuleDefinitions[mid].SetFirstDaqChannel(addedChannels);
            } else {
                fModuleDefinitions[mid].SetFirstDaqChannel(firstDaqChannel);
            }
 
            string decodingFile = GetFieldValue("decodingFile", moduleDefinition);
            fModuleDefinitions[mid].SetDecodingFile(decodingFile);
            addedChannels += fModuleDefinitions[mid].GetNumberOfChannels();
 
            moduleVector.push_back(std::move(fModuleDefinitions[mid]));
            // plane.AddModule( fModuleDefinitions[mid] );
 
            moduleDefinition = GetNextElement(moduleDefinition);
        }
 
        // We removed the constraint that module id's should start by 0 and have no
        // missing numbers in a multi-module readout plane. Modules can have their
        // special "id", e.g. M0, M2, M3, M4 in SJTU proto. We don't have M1
 
        for (auto& mod : moduleVector) {
            plane.AddModule(mod);
        }
 
        this->AddReadoutPlane(std::move(plane));
        planeDefinition = GetNextElement(planeDefinition);
    }
 
    ValidateReadout();
}
 
TRestDetectorReadoutModule* TRestDetectorReadout::ParseModuleDefinition(TiXmlElement* moduleDefinition) {
    auto mod = new TRestDetectorReadoutModule();
    TRestDetectorReadoutModule& module = *mod;
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Warning) module.EnableWarnings();
 
    module.SetName(GetFieldValue("name", moduleDefinition));
    module.SetSize(StringTo2DVector(GetFieldValue("size", moduleDefinition)));
    module.SetTolerance(StringToDouble(GetFieldValue("tolerance", moduleDefinition)));
    Double_t pixelTolerance = StringToDouble(GetFieldValue("pixelTolerance", moduleDefinition));
    if (pixelTolerance == -1) pixelTolerance = 1.e-6;
 
    vector<TRestDetectorReadoutChannel> channelVector;
    vector<int> channelIDVector;
    TiXmlElement* channelDefinition = GetElement("readoutChannel", moduleDefinition);
    while (channelDefinition != nullptr) {
        TRestDetectorReadoutChannel channel;
 
        Int_t id = StringToInteger(GetFieldValue("id", channelDefinition));
        if (id != -1) {
            channelIDVector.push_back(id);
        }
        channel.SetDaqID(-1);
 
        const string channelName = GetFieldValue("name", channelDefinition);
        channel.SetChannelName(channelName);
 
        vector<TRestDetectorReadoutPixel> pixelVector;
        vector<int> pixelIDVector;
        TiXmlElement* pixelDefinition = GetElement("addPixel", channelDefinition);
        while (pixelDefinition != nullptr) {
            TRestDetectorReadoutPixel pixel;
 
            pixel.SetOrigin(StringTo2DVector(GetFieldValue("origin", pixelDefinition)));
            pixel.SetSize(StringTo2DVector(GetFieldValue("size", pixelDefinition)));
            pixel.SetRotation(StringToDouble(GetFieldValue("rotation", pixelDefinition)));
            pixel.SetTriangle(StringToBool(GetFieldValue("triangle", pixelDefinition)));
            pixel.SetTolerance(pixelTolerance);
 
            if (StringToInteger(GetFieldValue("id", pixelDefinition)) != -1)
                pixelIDVector.push_back(StringToInteger(GetFieldValue("id", pixelDefinition)));
            pixelVector.push_back(pixel);
            // channel.AddPixel( pixel );
            pixelDefinition = GetNextElement(pixelDefinition);
        }
 
        if (!pixelIDVector.empty() && pixelIDVector.size() != pixelVector.size()) {
            RESTError
                << "pixel id definition may be wrong! It must be coherent and starts from 0. Check your "
                   "readout module definition!"
                << RESTendl;
            exit(0);
        }
 
        // Creating the vector fReadoutPixel in the channel with pixels added in the
        // order of their ID.
        for (Int_t i(0); i < (Int_t)pixelVector.size(); i++) {
            for (Int_t j(0); j < (Int_t)pixelVector.size(); j++) {
                if (pixelIDVector[j] == i) {
                    channel.AddPixel(pixelVector[j]);
                    break;
                }
            }
        }
 
        if (channel.GetNumberOfPixels() != (int)pixelVector.size()) {
            RESTError << "pixel id definition may be wrong! check your "
                         "readout module definition!"
                      << RESTendl;
            exit(0);
        }
 
        channelVector.push_back(channel);
        channelDefinition = GetNextElement(channelDefinition);
    }
 
    if (!channelIDVector.empty() && channelIDVector.size() != channelVector.size()) {
        RESTError << "TRestDetectorReadout::ParseModuleDefinition. Channel id definition may be wrong!"
                  << "check your readout module definition!" << RESTendl;
        RESTError << " " << RESTendl;
        RESTError << "channelIDVector size : " << channelIDVector.size() << RESTendl;
        RESTError << "channel vector size : " << channelVector.size() << RESTendl;
 
        exit(0);
    }
 
    // Creating the vector fReadoutChannel in the module with channels added in
    // the order of their ID.
    for (Int_t i(0); i < (Int_t)channelVector.size(); i++) {
        for (Int_t j(0); j < (Int_t)channelVector.size(); j++) {
            if (channelIDVector[j] == i) {
                module.AddChannel(channelVector[j]);
                break;
            }
        }
    }
 
    if (module.GetNumberOfChannels() != channelVector.size()) {
        RESTError << "TRestDetectorReadout::ParseModuleDefinition. Channel id definition may be wrong!"
                  << "check your readout module definition!" << RESTendl;
        RESTError << " " << RESTendl;
        RESTError << "Module number of channels : " << module.GetNumberOfChannels() << RESTendl;
        RESTError << "channel vector size : " << channelVector.size() << RESTendl;
 
        exit(0);
    }
 
    return mod;
}
 
void TRestDetectorReadout::ValidateReadout() const {
    RESTDebug << "--------------------------------------------------" << RESTendl;
    RESTDebug << "TRestDetectorReadout::ValidateReadout:: NOT IMPLEMENTED" << RESTendl;
    RESTDebug << "This function should crosscheck that there are no repeated "
                 "DaqChannels IDs"
              << RESTendl;
    RESTDebug << "If any checks are implemented in the future. Those checks should be "
                 "an option."
              << RESTendl;
    RESTDebug << "No dead area in the readout module" << RESTendl;
    RESTDebug << "And other checks" << RESTendl;
    RESTDebug << "--------------------------------------------------" << RESTendl;
}
 
void TRestDetectorReadout::GetPlaneModuleChannel(Int_t signalID, Int_t& planeID, Int_t& moduleID,
                                                 Int_t& channelID) {
    for (int p = 0; p < GetNumberOfReadoutPlanes(); p++) {
        TRestDetectorReadoutPlane* plane = &fReadoutPlanes[p];
        for (size_t m = 0; m < plane->GetNumberOfModules(); m++) {
            TRestDetectorReadoutModule* mod = &(*plane)[m];
 
            if (mod->IsDaqIDInside(signalID)) {
                planeID = plane->GetID();
                moduleID = mod->GetModuleID();
                channelID = mod->DaqToReadoutChannel(signalID);
            }
        }
    }
}
 
Int_t TRestDetectorReadout::GetHitsDaqChannel(const TVector3& position, Int_t& planeID, Int_t& moduleID,
                                              Int_t& channelID) {
    for (int p = 0; p < GetNumberOfReadoutPlanes(); p++) {
        TRestDetectorReadoutPlane* plane = &fReadoutPlanes[p];
        int m = plane->GetModuleIDFromPosition(position);
        if (m >= 0) {
            // TRestDetectorReadoutModule* mod = plane->GetModuleByID(m);
            TRestDetectorReadoutModule* mod = plane->GetModuleByID(m);
            Int_t readoutChannel = mod->FindChannel({position.X(), position.Y()});
            if (readoutChannel >= 0) {
                planeID = plane->GetID();
                moduleID = mod->GetModuleID();
                channelID = readoutChannel;
                return mod->GetChannel(readoutChannel)->GetDaqID();
            }
        }
    }
    return -1;
}
 
//
std::tuple<Int_t, Int_t, Int_t> TRestDetectorReadout::GetHitsDaqChannelAtReadoutPlane(
    const TVector3& position, Int_t planeId) {
    if (planeId > GetNumberOfReadoutPlanes()) {
        RESTWarning << "TRestDetectorReadout. Fail trying to retrieve planeId : " << planeId << RESTendl;
        RESTWarning << "Number of readout planes: " << GetNumberOfReadoutPlanes() << RESTendl;
        return std::make_tuple(-1, -1, -1);
    }
 
    set<Int_t> daqIds;
    Int_t moduleID = -1;
    Int_t channelID = -1;
 
    TRestDetectorReadoutPlane& plane = fReadoutPlanes[planeId];
    int m = plane.GetModuleIDFromPosition(position);
    if (m >= 0) {
        TRestDetectorReadoutModule* mod = plane.GetModuleByID(m);
 
        TRestDetectorReadoutChannel* channel = nullptr;
        int channelIndex = -1;
        if (mod->GetNumberOfChannels() == 1) {
            // workaround for vetoes which only have one channel
            channelIndex = 0;
            channel = mod->GetChannel(channelIndex);
        } else {
            channelIndex = mod->FindChannel({position.X(), position.Y()});
            channel = mod->GetChannel(channelIndex);
        }
        if (channel != nullptr) {
            moduleID = mod->GetModuleID();
            channelID = channelIndex;
            daqIds.insert(channel->GetDaqID());
        }
    }
 
    if (daqIds.empty()) {
        return std::make_tuple(-1, -1, -1);
    } else if (daqIds.size() == 1) {
        Int_t daqId = *daqIds.begin();
        return std::make_tuple(daqId, moduleID, channelID);
    } else {
        cerr << "TRestDetectorReadout::GetHitsDaqChannelAtReadoutPlane. More than one daq channel found for "
                "the given position. This means there is a problem with the readout definition."
             << endl;
        exit(1);
    }
}
 
Double_t TRestDetectorReadout::GetX(Int_t signalID) {
    Int_t planeID, readoutChannel = -1, readoutModule;
    GetPlaneModuleChannel(signalID, planeID, readoutModule, readoutChannel);
    if (readoutChannel == -1) {
        return std::numeric_limits<Double_t>::quiet_NaN();
    }
    return GetX(planeID, readoutModule, readoutChannel);
}
 
Double_t TRestDetectorReadout::GetY(Int_t signalID) {
    Int_t planeID, readoutChannel = -1, readoutModule;
    GetPlaneModuleChannel(signalID, planeID, readoutModule, readoutChannel);
    if (readoutChannel == -1) {
        return std::numeric_limits<Double_t>::quiet_NaN();
    }
    return GetY(planeID, readoutModule, readoutChannel);
}
 
Double_t TRestDetectorReadout::GetX(Int_t plane, Int_t modID, Int_t chID) {
    return GetReadoutPlaneWithID(plane)->GetX(modID, chID);
}
 
Double_t TRestDetectorReadout::GetY(Int_t plane, Int_t modID, Int_t chID) {
    return GetReadoutPlaneWithID(plane)->GetY(modID, chID);
}
 
void TRestDetectorReadout::PrintMetadata(Int_t DetailLevel) {
    if (DetailLevel >= 0) {
        TRestMetadata::PrintMetadata();
 
        RESTMetadata << "Number of readout planes : " << GetNumberOfReadoutPlanes() << RESTendl;
        RESTMetadata << "-----------------------------------" << RESTendl;
        for (int p = 0; p < GetNumberOfReadoutPlanes(); p++) {
            fReadoutPlanes[p].Print(DetailLevel - 1);
        }
        RESTMetadata << "****************************************" << RESTendl;
        cout << endl;
    }
}
 
void TRestDetectorReadout::Draw() {
    cout << " TRestDetectorReadout::Draw() is not implemented" << endl;
    cout << " To draw a TRestDetectorReadout class with name \"readoutName\"";
    cout << " stored in a ROOT file \"rootFile.root\"" << endl;
    cout << " You can use the script : REST_Readout_Viewer( \"rootFile.root\", "
            "\"readoutName\" )"
         << endl;
    cout << endl;
    cout << " Or you can access directly a readout plane and draw using : " << endl;
    cout << " readout->GetReadoutPlane( 0 )->Draw( ); " << endl;
}
 
void TRestDetectorReadout::Export(const string& fileName) {
    if (TRestTools::GetFileNameExtension(fileName) == "root") {
        TFile* f = TFile::Open(fileName.c_str(), "UPDATE");
        this->Write();
        f->Close();
    } else {
        RESTWarning << "Can only export readout as a root file, skipping..." << RESTendl;
    }
}
 
Int_t TRestDetectorReadout::GetDaqId(const TVector3& position, bool check) {
    std::vector<int> daqIds;
    for (int planeIndex = 0; planeIndex < GetNumberOfReadoutPlanes(); planeIndex++) {
        // const TRestDetectorReadoutPlane& plane = fReadoutPlanes[planeIndex];
        const auto [daqId, moduleID, channelID] = GetHitsDaqChannelAtReadoutPlane(position, planeIndex);
        if (daqId != -1) {
            if (!check) {
                return daqId;
            }
            daqIds.push_back(daqId);
        }
    }
 
    if (daqIds.empty()) {
        return -1;
    } else if (daqIds.size() == 1) {
        return daqIds[0];
    } else {
        cerr << "TRestDetectorReadout::GetDaqId. More than one daq channel found for "
                "the given position. This means there is a problem with the readout definition."
             << endl;
        exit(1);
    }
}
 
set<Int_t> TRestDetectorReadout::GetAllDaqIds() {
    set<Int_t> daqIds;
 
    for (int p = 0; p < GetNumberOfReadoutPlanes(); p++) {
        for (size_t m = 0; m < fReadoutPlanes[p].GetNumberOfModules(); m++) {
            for (size_t c = 0; c < fReadoutPlanes[p][m].GetNumberOfChannels(); c++) {
                auto channel = fReadoutPlanes[p][m][c];
                daqIds.insert(channel.GetDaqID());
            }
        }
    }
 
    return daqIds;
}
 
string TRestDetectorReadout::GetTypeForChannelDaqId(Int_t daqId) {
    for (int p = 0; p < GetNumberOfReadoutPlanes(); p++) {
        for (size_t m = 0; m < fReadoutPlanes[p].GetNumberOfModules(); m++) {
            for (size_t c = 0; c < fReadoutPlanes[p][m].GetNumberOfChannels(); c++) {
                auto channel = fReadoutPlanes[p][m][c];
                if (channel.GetDaqID() == daqId) {
                    return channel.GetType();
                }
            }
        }
    }
    return {};
}