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"
 
bool RESTREADOUT_DECODINGFILE_ERROR = false;
 
ClassImp(TRestDetectorReadout);
TRestDetectorReadout::TRestDetectorReadout() { Initialize(); }
 
TRestDetectorReadout::TRestDetectorReadout(const char* configFilename) : TRestMetadata(configFilename) {
    std::cout << "Loading readout. This might take few seconds" << std::endl;
    Initialize();
 
    LoadConfigFromFile(fConfigFileName);
}
 
TRestDetectorReadout::TRestDetectorReadout(const char* configFilename, std::string name)
    : TRestMetadata(configFilename) {
    std::cout << "Loading readout. This might take few seconds" << std::endl;
    Initialize();
 
    LoadConfigFromFile(fConfigFileName, name);
}
 
void TRestDetectorReadout::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fDecoding = false;
 
    fNReadoutPlanes = 0;
    fReadoutPlanes.clear();
}
 
TRestDetectorReadout::~TRestDetectorReadout() {}
 
Int_t TRestDetectorReadout::GetNumberOfReadoutPlanes() { return fReadoutPlanes.size(); }
 
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(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) {
    int planeID = -1, moduleID = -1, channelID = -1;
 
    GetPlaneModuleChannel(daqId, planeID, moduleID, channelID);
    if (channelID >= 0) {
        return &fReadoutPlanes[planeID][moduleID][channelID];
    }
 
    return nullptr;
}
 
TRestDetectorReadoutPlane* TRestDetectorReadout::GetReadoutPlane(int p) {
    if (p < fNReadoutPlanes)
        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 : " << fNReadoutPlanes << RESTendl;
    }
 
    return nullptr;
}
 
void TRestDetectorReadout::AddReadoutPlane(TRestDetectorReadoutPlane plane) {
    fReadoutPlanes.push_back(plane);
    fNReadoutPlanes++;
}
 
void TRestDetectorReadout::InitFromConfigFile() {
    fMappingNodes = StringToInteger(GetParameter("mappingNodes", "0"));
 
    TiXmlElement* moduleDefinition = GetElement("readoutModule");
    while (moduleDefinition != nullptr) {
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
            std::cout << "------module-----------------" << std::endl;
            std::cout << moduleDefinition << std::endl;
            std::cout << "-----------------------------" << std::endl;
            GetChar();
        }
 
        TRestDetectorReadoutModule module = *ParseModuleDefinition(moduleDefinition);
        module.DoReadoutMapping(fMappingNodes);
        fModuleDefinitions.push_back(module);
        moduleDefinition = GetNextElement(moduleDefinition);
    }
 
    TiXmlElement* planeDefinition = GetElement("readoutPlane");
    std::vector<TRestDetectorReadoutModule> moduleVector;
    Int_t addedChannels = 0;
    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));
 
        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) firstDaqChannel = addedChannels;
 
            std::string decodingFile = GetFieldValue("decodingFile", moduleDefinition);
            if (decodingFile == "Not defined" || decodingFile.empty() || RESTREADOUT_DECODINGFILE_ERROR)
                fDecoding = false;
            else
                fDecoding = true;
 
            if (fDecoding && !TRestTools::fileExists(decodingFile)) {
                RESTWarning << "The decoding file does not exist!" << RESTendl;
                RESTWarning << "--------------------------------" << RESTendl;
                RESTWarning << "File : " << decodingFile << RESTendl;
                RESTWarning << "Default decoding will be used. readoutChannel=daqChannel" << RESTendl;
                RESTWarning << "To avoid this message and use the default decoding define : "
                               "decodingFile=\"\""
                            << RESTendl;
                RESTWarning << "--------------------------------" << RESTendl;
                RESTWarning << "Press a KEY to continue..." << RESTendl;
                getchar();
                fDecoding = false;
                RESTREADOUT_DECODINGFILE_ERROR = true;
            }
 
            std::vector<Int_t> rChannel;
            std::vector<Int_t> dChannel;
            if (fDecoding && TRestTools::fileExists(decodingFile)) {
                FILE* f = fopen(decodingFile.c_str(), "r");
 
                Int_t daq, readout;
                while (!feof(f)) {
                    if (fscanf(f, "%d\t%d\n", &daq, &readout) <= 0) {
                        RESTError << "TRestDetectorReadout::InitFromConfigFile. Problem reading decoding"
                                  << RESTendl;
                        RESTError << "This error might need support at REST forum" << RESTendl;
                        exit(-1);
                    }
                    // we skip blank daq channels if readout id is <0
                    // e.g. daq id: 22, readout id: -1
                    if (readout >= 0) {
                        rChannel.push_back(readout);
                        dChannel.push_back(daq + firstDaqChannel);
                    }
                }
                fclose(f);
            }
 
            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
                std::cout << "------module-----------------" << std::endl;
                std::cout << moduleDefinition << std::endl;
                std::cout << "-----------------------------" << std::endl;
                getchar();
            }
 
            if (fDecoding && (unsigned int)fModuleDefinitions[mid].GetNumberOfChannels() != rChannel.size()) {
                RESTError << "TRestDetectorReadout."
                          << " The number of channels defined in the readout is not the same"
                          << " as the number of channels found in the decoding." << RESTendl;
                exit(1);
            }
 
            for (size_t ch = 0; ch < fModuleDefinitions[mid].GetNumberOfChannels(); ch++) {
                if (!fDecoding) {
                    Int_t id = ch;
                    rChannel.push_back(id);
                    dChannel.push_back(id + firstDaqChannel);
                }
 
                // WRONG version before -->
                // fModuleDefinitions[mid].GetChannel(ch)->SetID( rChannel[ch] );
                if (!fModuleDefinitions[mid].GetChannel(rChannel[ch])) {
                    RESTError << "Problem setting readout channel " << rChannel[ch]
                              << " with daq id: " << dChannel[ch] << RESTendl;
                    continue;
                }
                fModuleDefinitions[mid].GetChannel(rChannel[ch])->SetDaqID(dChannel[ch]);
                fModuleDefinitions[mid].GetChannel(rChannel[ch])->SetChannelID(rChannel[ch]);
 
                addedChannels++;
            }
            fModuleDefinitions[mid].SetMinMaxDaqIDs();
 
            moduleVector.push_back(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& i : moduleVector) {
            plane.AddModule(i);
        }
 
        this->AddReadoutPlane(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;
 
    std::vector<TRestDetectorReadoutChannel> channelVector;
    std::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);
 
        std::vector<TRestDetectorReadoutPixel> pixelVector;
        std::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 std::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 std::vector size : " << channelVector.size() << RESTendl;
 
        exit(0);
    }
 
    // Creating the std::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 std::vector size : " << channelVector.size() << RESTendl;
 
        exit(0);
    }
 
    return mod;
}
 
void TRestDetectorReadout::ValidateReadout() {
    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.X(), position.Y(), position.Z());
        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;
}
 
//
Int_t TRestDetectorReadout::GetHitsDaqChannelAtReadoutPlane(const TVector3& hitpos, Int_t& moduleID,
                                                            Int_t& channelID, Int_t planeId) {
    if (planeId > GetNumberOfReadoutPlanes()) {
        RESTWarning << "TRestDetectorReadout. Fail trying to retrieve planeId : " << planeId << RESTendl;
        RESTWarning << "Number of readout planes: " << GetNumberOfReadoutPlanes() << RESTendl;
        return -1;
    }
 
    TRestDetectorReadoutPlane* plane = &fReadoutPlanes[planeId];
    int m = plane->GetModuleIDFromPosition(hitpos.X(), hitpos.Y(), hitpos.Z());
    if (m >= 0) {
        TRestDetectorReadoutModule* mod = plane->GetModuleByID(m);
        Int_t readoutChannel = mod->FindChannel({hitpos.X(), hitpos.Y()});
        if (readoutChannel >= 0) {
            moduleID = mod->GetModuleID();
            channelID = readoutChannel;
            return mod->GetChannel(readoutChannel)->GetDaqID();
        }
    }
    return -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 : " << fNReadoutPlanes << RESTendl;
        RESTMetadata << "Decoding was defined : ";
        if (fDecoding)
            RESTMetadata << "YES" << RESTendl;
        else
            RESTMetadata << "NO" << RESTendl;
        RESTMetadata << "-----------------------------------" << RESTendl;
        for (int p = 0; p < GetNumberOfReadoutPlanes(); p++) fReadoutPlanes[p].Print(DetailLevel - 1);
        RESTMetadata << "****************************************" << RESTendl;
        std::cout << std::endl;
    }
}
 
void TRestDetectorReadout::Draw() {
    std::cout << " TRestDetectorReadout::Draw() is not implemented" << std::endl;
    std::cout << " To draw a TRestDetectorReadout class with name \"readoutName\"";
    std::cout << " stored in a ROOT file \"rootFile.root\"" << std::endl;
    std::cout << " You can use the script : REST_Readout_Viewer( \"rootFile.root\", "
                 "\"readoutName\" )"
              << std::endl;
    std::cout << std::endl;
    std::cout << " Or you can access directly a readout plane and draw using : " << std::endl;
    std::cout << " readout->GetReadoutPlane( 0 )->Draw( ); " << std::endl;
}