TRestDetectorHitsGaussAnalysisProcess.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 "TRestDetectorHitsGaussAnalysisProcess.h"
 
using namespace std;
 
ClassImp(TRestDetectorHitsGaussAnalysisProcess);
 
TRestDetectorHitsGaussAnalysisProcess::TRestDetectorHitsGaussAnalysisProcess() { Initialize(); }
 
TRestDetectorHitsGaussAnalysisProcess::TRestDetectorHitsGaussAnalysisProcess(const char* configFilename) {
    Initialize();
 
    LoadConfig(configFilename);
}
 
TRestDetectorHitsGaussAnalysisProcess::~TRestDetectorHitsGaussAnalysisProcess() {
    delete fInputHitsEvent;
    delete fOutputHitsEvent;
}
 
void TRestDetectorHitsGaussAnalysisProcess::LoadDefaultConfig() {
    SetName(this->ClassName());
    SetTitle("Default config");
}
 
void TRestDetectorHitsGaussAnalysisProcess::LoadConfig(const string& configFilename, const string& name) {
    if (LoadConfigFromFile(configFilename, name)) LoadDefaultConfig();
}
 
void TRestDetectorHitsGaussAnalysisProcess::InitProcess() {
    // For example, try to initialize a pointer to existing metadata
    // accessible from TRestRun
}
 
void TRestDetectorHitsGaussAnalysisProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    // fHitsEvent = new TRestDetectorHitsEvent();
    fInputHitsEvent = nullptr;
    fOutputHitsEvent = new TRestDetectorHitsEvent();
}
 
TRestEvent* TRestDetectorHitsGaussAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) {
    fInputHitsEvent = (TRestDetectorHitsEvent*)inputEvent;
 
    TString obsName;
 
    // cout << "Event id : " << fInputHitsEvent->GetID() << endl;
    TRestHits* hits = fInputHitsEvent->GetHits();
 
    for (unsigned int n = 0; n < hits->GetNumberOfHits(); n++) {
        Double_t eDep = hits->GetEnergy(n);
 
        Double_t x = hits->GetX(n);
        Double_t y = hits->GetY(n);
        Double_t z = hits->GetZ(n);
 
        auto time = hits->GetTime(n);
        auto type = hits->GetType(n);
 
        fOutputHitsEvent->AddHit(x, y, z, eDep, time, type);
    }
 
    Double_t gausSigmaX = fOutputHitsEvent->GetGaussSigmaX(fError, fNHitsMin);
    Double_t gausSigmaY = fOutputHitsEvent->GetGaussSigmaY(fError, fNHitsMin);
    Double_t gausSigmaZ = fOutputHitsEvent->GetGaussSigmaZ(fError, fNHitsMin);
    Double_t xy2SigmaGaus = (gausSigmaX == -1. || gausSigmaY == -1.)
                                ? -1.
                                : (gausSigmaX * gausSigmaX) + (gausSigmaY * gausSigmaY);
    Double_t xySigmaBalanceGaus = (gausSigmaX == -1. || gausSigmaY == -1.)
                                      ? -1.
                                      : (gausSigmaX - gausSigmaY) / (gausSigmaX + gausSigmaY);
 
    if (hits->GetNumberOfHits() > 30 && xy2SigmaGaus < 0.05)
        RESTDebug << string("Event ID: ") << to_string(fInputHitsEvent->GetID()) << string("||") << RESTendl;
 
    SetObservableValue("xSigmaGaus", gausSigmaX);
    SetObservableValue("ySigmaGaus", gausSigmaY);
    SetObservableValue("zSigmaGaus", gausSigmaZ);
    SetObservableValue("xy2SigmaGaus", xy2SigmaGaus);
    SetObservableValue("xySigmaBalanceGaus", xySigmaBalanceGaus);
 
    // We transform here fHitsEvent if necessary
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        fOutputHitsEvent->PrintEvent();
 
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) GetChar();
    }
 
    return fOutputHitsEvent;
}
 
void TRestDetectorHitsGaussAnalysisProcess::PrintMetadata() {
    BeginPrintProcess();
 
    // Print output metadata using, metadata << endl;
    RESTMetadata << "Pitch (mm) : " << fPitch << RESTendl;
    RESTMetadata << "Error (ADC) : " << fError << RESTendl;
    RESTMetadata << "Minimum number of hits to apply correction : " << fNHitsMin << RESTendl;
 
    EndPrintProcess();
}