TRestDetectorHitsToSignalProcess.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 "TRestDetectorHitsToSignalProcess.h"
 
using namespace std;
 
ClassImp(TRestDetectorHitsToSignalProcess);
 
TRestDetectorHitsToSignalProcess::TRestDetectorHitsToSignalProcess() { Initialize(); }
 
TRestDetectorHitsToSignalProcess::TRestDetectorHitsToSignalProcess(const char* configFilename) {
    Initialize();
 
    if (LoadConfigFromFile(configFilename) == -1) {
        LoadDefaultConfig();
    }
 
    if (fReadout == nullptr) {
        fReadout = new TRestDetectorReadout(configFilename);
    }
}
 
TRestDetectorHitsToSignalProcess::~TRestDetectorHitsToSignalProcess() {
    delete fReadout;
    delete fSignalEvent;
}
 
void TRestDetectorHitsToSignalProcess::LoadDefaultConfig() {
    SetName("hitsToSignalProcess-Default");
    SetTitle("Default config");
 
    cout << "Hits to signal metadata not found. Loading default values" << endl;
 
    fSampling = 1;
    fElectricField = -1;
    fGasPressure = -1;
}
 
void TRestDetectorHitsToSignalProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fReadout = nullptr;
    fGas = nullptr;
 
    fHitsEvent = nullptr;
    fSignalEvent = new TRestDetectorSignalEvent();
}
 
void TRestDetectorHitsToSignalProcess::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 "
               "TRestDetectorHitsToSignalProcess"
            << 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();
        }
    }
 
    if (fDriftVelocity <= 0) {
        RESTError
            << "Drift velocity not defined. Please, define it in the process metadata or in TRestDetectorGas"
            << RESTendl;
        exit(1);
    }
 
    fReadout = GetMetadata<TRestDetectorReadout>();
    if (fReadout == nullptr) {
        if (!this->GetError()) {
            this->SetError("The readout was not properly initialized.");
        }
    }
}
 
TRestEvent* TRestDetectorHitsToSignalProcess::ProcessEvent(TRestEvent* inputEvent) {
    fHitsEvent = (TRestDetectorHitsEvent*)inputEvent;
    fSignalEvent->SetEventInfo(fHitsEvent);
 
    if (!fReadout) {
        return nullptr;
    }
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        cout << "Number of hits : " << fHitsEvent->GetNumberOfHits() << endl;
        cout << "--------------------------" << endl;
    }
 
    for (unsigned int hit = 0; hit < fHitsEvent->GetNumberOfHits(); hit++) {
        Double_t x = fHitsEvent->GetX(hit);
        Double_t y = fHitsEvent->GetY(hit);
        Double_t z = fHitsEvent->GetZ(hit);
        Double_t t = fHitsEvent->GetTime(hit);
 
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme && hit < 20) {
            cout << "Hit : " << hit << " x : " << x << " y : " << y << " z : " << z << " t : " << t << endl;
        }
 
        for (int p = 0; p < fReadout->GetNumberOfReadoutPlanes(); p++) {
            const auto [daqId, moduleId, channelId] = fReadout->GetHitsDaqChannelAtReadoutPlane({x, y, z}, p);
            TRestDetectorReadoutPlane* plane = fReadout->GetReadoutPlane(p);
 
            if (daqId >= 0) {
                auto channel = fReadout->GetReadoutChannelWithDaqID(daqId);
 
                const bool isVeto = (channel->GetChannelType() == "veto");
 
                Double_t energy = fHitsEvent->GetEnergy(hit);
                const auto distance = plane->GetDistanceTo({x, y, z});
                if (distance < 0) {
                    RESTError << "TRestDetectorHitsToSignalProcess: Negative distance to readout plane. "
                                 "This should not happen."
                              << RESTendl;
                    exit(1);
                }
 
                auto velocity = fDriftVelocity;
                if (isVeto) {
                    velocity = REST_Physics::lightSpeed;
                }
 
                if (velocity <= 0) {
                    RESTError
                        << "TRestDetectorHitsToSignalProcess: Negative velocity. This should not happen."
                        << RESTendl;
                    exit(1);
                }
 
                Double_t time = t + distance / velocity;
 
                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug && hit < 20) {
                    cout << "Module : " << moduleId << " Channel : " << channelId << " daq ID : " << daqId
                         << endl;
                }
 
                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug && hit < 20) {
                    cout << "Energy : " << energy << " time : " << time << endl;
                }
                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme && hit < 20) {
                    printf(
                        " TRestDetectorHitsToSignalProcess: x %lf y %lf z %lf energy %lf t %lf "
                        "fDriftVelocity %lf fSampling %lf time %lf\n",
                        x, y, z, energy, t, fDriftVelocity, fSampling, time);
                }
 
                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
                    cout << "Drift velocity : " << fDriftVelocity << " mm/us" << endl;
                }
 
                time = floor(time / fSampling) * fSampling;
 
                if (time < 0) {
                    RESTError << "TRestDetectorHitsToSignalProcess: Negative time. This should not happen. "
                                 "EventID: "
                              << fHitsEvent->GetID() << RESTendl;
                    exit(1);
                }
                fSignalEvent->AddChargeToSignal(daqId, time, energy);
 
                auto signal = fSignalEvent->GetSignalById(daqId);
                signal->SetSignalName(channel->GetChannelName());
                signal->SetSignalType(channel->GetChannelType());
 
            } else {
                if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
                    RESTDebug << "TRestDetectorHitsToSignalProcess. Readout channel not find for position ("
                              << x << ", " << y << ", " << z << ")!" << RESTendl;
                }
            }
        }
    }
 
    fSignalEvent->SortSignals();
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        cout << "TRestDetectorHitsToSignalProcess : Number of signals added : "
             << fSignalEvent->GetNumberOfSignals() << endl;
        cout << "TRestDetectorHitsToSignalProcess : Total signals integral : " << fSignalEvent->GetIntegral()
             << endl;
    }
 
    if (fSignalEvent->GetNumberOfSignals() == 0) {
        return nullptr;
    }
 
    return fSignalEvent;
}