TRestDetectorSignalToHitsProcess.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 "TRestDetectorSignalToHitsProcess.h"
 
#include "TRestDetectorSetup.h"
 
using namespace std;
 
ClassImp(TRestDetectorSignalToHitsProcess);
 
TRestDetectorSignalToHitsProcess::TRestDetectorSignalToHitsProcess() { Initialize(); }
 
TRestDetectorSignalToHitsProcess::TRestDetectorSignalToHitsProcess(const char* configFilename) {
    Initialize();
 
    if (LoadConfigFromFile(configFilename) == -1) {
        LoadDefaultConfig();
    }
}
 
void TRestDetectorSignalToHitsProcess::LoadDefaultConfig() { SetTitle("Default config"); }
 
TRestDetectorSignalToHitsProcess::~TRestDetectorSignalToHitsProcess() { delete fHitsEvent; }
 
void TRestDetectorSignalToHitsProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fHitsEvent = new TRestDetectorHitsEvent();
    fSignalEvent = nullptr;
 
    fGas = nullptr;
    fReadout = nullptr;
}
 
void TRestDetectorSignalToHitsProcess::InitProcess() {
    fGas = GetMetadata<TRestDetectorGas>();
    if (fGas != nullptr) {
#ifndef USE_Garfield
        RESTError << "A TRestDetectorGas definition was found but REST was not linked to Garfield libraries."
                  << RESTendl;
        RESTError
            << "Please, remove the TRestDetectorGas definition, and add gas parameters inside the process "
               "TRestDetectorSignalToHitsProcess"
            << RESTendl;
        if (!fGas->GetError()) fGas->SetError("REST was not compiled with Garfield.");
        if (!this->GetError()) this->SetError("Attempt to use TRestDetectorGas without Garfield");
#endif
 
        if (fGasPressure <= 0) {
            fGasPressure = fGas->GetPressure();
        }
        if (fElectricField <= 0) {
            fElectricField = fGas->GetElectricField();
        }
 
        fGas->SetPressure(fGasPressure);
        fGas->SetElectricField(fElectricField);
 
        if (fDriftVelocity <= 0) {
            fDriftVelocity = fGas->GetDriftVelocity();
        }
    } else {
        if (fDriftVelocity < 0) {
            if (!this->GetError()) {
                this->SetError("Drift velocity is negative.");
            }
        }
    }
 
    fReadout = GetMetadata<TRestDetectorReadout>();
 
    if (fReadout == nullptr) {
        if (!this->GetError()) {
            this->SetError("The readout was not properly initialized.");
        }
    }
}
 
TRestEvent* TRestDetectorSignalToHitsProcess::ProcessEvent(TRestEvent* inputEvent) {
    fSignalEvent = (TRestDetectorSignalEvent*)inputEvent;
 
    if (!fReadout) {
        return nullptr;
    }
 
    fHitsEvent->SetID(fSignalEvent->GetID());
    fHitsEvent->SetSubID(fSignalEvent->GetSubID());
    fHitsEvent->SetTimeStamp(fSignalEvent->GetTimeStamp());
    fHitsEvent->SetSubEventTag(fSignalEvent->GetSubEventTag());
 
    RESTDebug << "TRestDetectorSignalToHitsProcess. Event id : " << fHitsEvent->GetID() << RESTendl;
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
        fSignalEvent->PrintEvent();
    }
 
    Int_t numberOfSignals = fSignalEvent->GetNumberOfSignals();
 
    if (numberOfSignals == 0) return nullptr;
 
    Int_t planeID, readoutChannel = -1, readoutModule;
    for (int i = 0; i < numberOfSignals; i++) {
        TRestDetectorSignal* signal = fSignalEvent->GetSignal(i);
        Int_t signalID = signal->GetSignalID();
 
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug)
            cout << "Searching readout coordinates for signal ID : " << signalID << endl;
 
        fReadout->GetPlaneModuleChannel(signalID, planeID, readoutModule, readoutChannel);
 
        if (readoutChannel == -1) {
            // cout << "REST Warning : Readout channel not found for daq ID : " << signalID << endl;
            continue;
        }
 
        TRestDetectorReadoutPlane* plane = fReadout->GetReadoutPlaneWithID(planeID);
 
        // For the moment this will only be valid for a TPC with its axis (field
        // direction) being in z
        Double_t fieldZDirection = plane->GetNormal().Z();
        Double_t zPosition = plane->GetPosition().Z();
 
        Double_t x = plane->GetX(readoutModule, readoutChannel);
        Double_t y = plane->GetY(readoutModule, readoutChannel);
 
        REST_HitType type = XYZ;
        TRestDetectorReadoutModule* mod = plane->GetModuleByID(readoutModule);
        if (TMath::IsNaN(x)) {
            x = mod->GetPlaneCoordinates(TVector2(mod->GetSize().X() / 2, mod->GetSize().Y() / 2)).X();
            type = YZ;
        } else if (TMath::IsNaN(y)) {
            y = mod->GetPlaneCoordinates(TVector2(mod->GetSize().X() / 2, mod->GetSize().Y() / 2)).Y();
            type = XZ;
        }
 
        if (fMethod == "onlyMax") {
            Double_t hitTime = signal->GetMaxPeakTime();
            Double_t distanceToPlane = hitTime * fDriftVelocity;
 
            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug)
                cout << "Distance to plane : " << distanceToPlane << endl;
 
            Double_t z = zPosition + fieldZDirection * distanceToPlane;
 
            Double_t energy = signal->GetMaxPeakValue();
 
            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug)
                cout << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << z
                     << " Energy : " << energy << endl;
 
            fHitsEvent->AddHit(x, y, z, energy, 0, type);
        } else if (fMethod == "tripleMax") {
            Int_t bin = signal->GetMaxIndex();
            int binprev = (bin - 1) < 0 ? bin : bin - 1;
            int binnext = (bin + 1) > signal->GetNumberOfPoints() - 1 ? bin : bin + 1;
 
            Double_t hitTime = signal->GetTime(bin);
            Double_t energy = signal->GetData(bin);
 
            Double_t distanceToPlane = hitTime * fDriftVelocity;
            Double_t z = zPosition + fieldZDirection * distanceToPlane;
 
            fHitsEvent->AddHit(x, y, z, energy, 0, type);
 
            hitTime = signal->GetTime(binprev);
            energy = signal->GetData(binprev);
 
            distanceToPlane = hitTime * fDriftVelocity;
            z = zPosition + fieldZDirection * distanceToPlane;
 
            fHitsEvent->AddHit(x, y, z, energy, 0, type);
 
            hitTime = signal->GetTime(binnext);
            energy = signal->GetData(binnext);
 
            distanceToPlane = hitTime * fDriftVelocity;
            z = zPosition + fieldZDirection * distanceToPlane;
 
            fHitsEvent->AddHit(x, y, z, energy, 0, type);
 
            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
                cout << "Distance to plane : " << distanceToPlane << endl;
                cout << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << z
                     << " Energy : " << energy << endl;
            }
        } else if (fMethod == "tripleMaxAverage") {
            Int_t bin = signal->GetMaxIndex();
            int binprev = (bin - 1) < 0 ? bin : bin - 1;
            int binnext = (bin + 1) > signal->GetNumberOfPoints() - 1 ? bin : bin + 1;
 
            Double_t hitTime = signal->GetTime(bin);
            Double_t energy1 = signal->GetData(bin);
 
            Double_t distanceToPlane = hitTime * fDriftVelocity;
            Double_t z1 = zPosition + fieldZDirection * distanceToPlane;
 
            hitTime = signal->GetTime(binprev);
            Double_t energy2 = signal->GetData(binprev);
 
            distanceToPlane = hitTime * fDriftVelocity;
            Double_t z2 = zPosition + fieldZDirection * distanceToPlane;
 
            hitTime = signal->GetTime(binnext);
            Double_t energy3 = signal->GetData(binnext);
 
            distanceToPlane = hitTime * fDriftVelocity;
            Double_t z3 = zPosition + fieldZDirection * distanceToPlane;
 
            Double_t eTot = energy1 + energy2 + energy3;
 
            Double_t zAvg = ((z1 * energy1) + (z2 * energy2) + (z3 * energy3)) / eTot;
            Double_t eAvg = eTot / 3.0;
 
            fHitsEvent->AddHit(x, y, zAvg, eAvg, 0, type);
 
            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
                cout << "Distance to plane: " << distanceToPlane << endl;
                cout << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << zAvg
                     << " Energy : " << eAvg << endl;
                cout << "z1, z2, z3 = " << z1 << ", " << z2 << ", " << z3 << endl;
                cout << "E1, E2, E3 = " << energy1 << ", " << energy2 << ", " << energy3 << endl;
            }
 
        } else if (fMethod == "gaussFit" || fMethod == "landauFit" || fMethod == "agetFit") {
            optional<pair<Double_t, Double_t>> peak;
            if (fMethod == "gaussFit") {
                peak = signal->GetPeakGauss();
            } else if (fMethod == "landauFit") {
                peak = signal->GetPeakLandau();
            } else if (fMethod == "agetFit") {
                peak = signal->GetPeakAget();
            } else {
                throw std::runtime_error("Invalid method");
            }
            if (!peak) {
                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) {
                    cout << "Unable to find peak for signal " << signal->GetSignalID()
                         << " with method: " << fMethod << endl;
                }
                continue;
            }
            const auto [time, energy] = peak.value();
 
            const auto distanceToPlane = time * fDriftVelocity;
            const auto z = zPosition + fieldZDirection * distanceToPlane;
 
            if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
                cout << "Signal event : " << signal->GetSignalID()
                     << "--------------------------------------------------------" << endl;
                cout << "Method: " << fMethod << " : time " << time << " us and energy : " << energy << endl;
                cout << "Signal to hit info : zPosition : " << zPosition
                     << "; fieldZDirection : " << fieldZDirection << " and driftV : " << fDriftVelocity
                     << endl;
                cout << "Adding hit. Time : " << time << " us x : " << x << " y : " << y << " z : " << z
                     << " Energy : " << energy << endl;
            }
 
            fHitsEvent->AddHit(x, y, z, energy, 0, type);
        } else if (fMethod == "qCenter") {
            Double_t energy_signal = 0;
            Double_t distanceToPlane = 0;
 
            for (int j = 0; j < signal->GetNumberOfPoints(); j++) {
                Double_t energy_point = signal->GetData(j);
                energy_signal += energy_point;
                distanceToPlane += signal->GetTime(j) * fDriftVelocity * energy_point;
            }
            Double_t energy = energy_signal / signal->GetNumberOfPoints();
 
            Double_t z = zPosition + fieldZDirection * (distanceToPlane / energy_signal);
            fHitsEvent->AddHit(x, y, z, energy, 0, type);
        } else if (fMethod == "all") {
            for (int j = 0; j < signal->GetNumberOfPoints(); j++) {
                Double_t energy = signal->GetData(j);
 
                Double_t distanceToPlane = signal->GetTime(j) * fDriftVelocity;
 
                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
                    cout << "Time : " << signal->GetTime(j) << " Drift velocity : " << fDriftVelocity << endl;
                    cout << "Distance to plane : " << distanceToPlane << endl;
                }
 
                Double_t z = zPosition + fieldZDirection * distanceToPlane;
 
                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug)
                    cout << "Adding hit. Time : " << signal->GetTime(j) << " x : " << x << " y : " << y
                         << " z : " << z << endl;
 
                fHitsEvent->AddHit(x, y, z, energy, 0, type);
            }
        } else if (fMethod == "intwindow") {
            Int_t nPoints = signal->GetNumberOfPoints();
            std::map<int, std::pair<int, double>> windowMap;
            for (int j = 0; j < nPoints; j++) {
                int index = signal->GetTime(j) / fIntWindow;
                auto it = windowMap.find(index);
                if (it != windowMap.end()) {
                    it->second.first++;
                    it->second.second += signal->GetData(j);
                } else {
                    windowMap[index] = std::make_pair(1, signal->GetData(j));
                }
            }
 
            for (const auto& [index, pair] : windowMap) {
                Double_t hitTime = index * fIntWindow + fIntWindow / 2.;
                Double_t energy = pair.second / pair.first;
                if (energy < fThreshold) {
                    continue;
                }
                RESTDebug << "TimeBin " << index << " Time " << hitTime << " Charge: " << energy
                          << " Thr: " << (fThreshold) << RESTendl;
                Double_t distanceToPlane = hitTime * fDriftVelocity;
                Double_t z = zPosition + fieldZDirection * distanceToPlane;
 
                RESTDebug << "Time : " << hitTime << " Drift velocity : " << fDriftVelocity
                          << "\nDistance to plane : " << distanceToPlane << RESTendl;
                RESTDebug << "Adding hit. Time : " << hitTime << " x : " << x << " y : " << y << " z : " << z
                          << " type " << type << RESTendl;
 
                fHitsEvent->AddHit(x, y, z, energy, 0, type);
            }
        } else {
            string errMsg = "The method " + (string)fMethod + " is not implemented!";
            SetError(errMsg);
        }
    }
 
    RESTDebug << "TRestDetectorSignalToHitsProcess. Hits added : " << fHitsEvent->GetNumberOfHits()
              << RESTendl;
    RESTDebug << "TRestDetectorSignalToHitsProcess. Hits total energy : " << fHitsEvent->GetTotalEnergy()
              << RESTendl;
 
    if (this->GetVerboseLevel() == TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        fHitsEvent->PrintEvent(30);
    } else if (this->GetVerboseLevel() == TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
        fHitsEvent->PrintEvent(-1);
    }
 
    if (fHitsEvent->GetNumberOfHits() <= 0) {
        string errMsg = "Last event id: " + IntegerToString(fHitsEvent->GetID()) +
                        ". Failed to find readout positions in channel to hit conversion.";
        SetWarning(errMsg);
        return nullptr;
    }
 
    return fHitsEvent;
}