TRestAxionFieldPropagationProcess.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 "TRestAxionFieldPropagationProcess.h"
 
#include <TVectorD.h>
 
#include <numeric>
using namespace std;
 
ClassImp(TRestAxionFieldPropagationProcess);
 
TRestAxionFieldPropagationProcess::TRestAxionFieldPropagationProcess() { Initialize(); }
 
TRestAxionFieldPropagationProcess::TRestAxionFieldPropagationProcess(char* cfgFileName) { Initialize(); }
 
TRestAxionFieldPropagationProcess::~TRestAxionFieldPropagationProcess() { delete fAxionEvent; }
 
void TRestAxionFieldPropagationProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fAxionEvent = new TRestAxionEvent();
}
 
void TRestAxionFieldPropagationProcess::InitProcess() {
    RESTDebug << "Entering ... TRestAxionFieldPropagationProcess::InitProcess" << RESTendl;
 
    fMagneticField = (TRestAxionMagneticField*)this->GetMetadata("TRestAxionMagneticField");
 
    if (!fMagneticField) {
        RESTError << "TRestAxionFieldPropagationprocess. Magnetic Field was not defined!" << RESTendl;
        exit(0);
    }
 
    for (size_t n = 0; n < fMagneticField->GetNumberOfVolumes(); n++) fMagneticField->ReMap(n, fReMap);
 
    if (!fAxionField) {
        fAxionField = new TRestAxionField();
 
        fBufferGas = (TRestAxionBufferGas*)this->GetMetadata("TRestAxionBufferGas");
        if (fBufferGas)
            fAxionField->AssignBufferGas(fBufferGas);
        else
            fBufferGasAdditionalLength = 0;
 
        fAxionField->AssignMagneticField(fMagneticField);
        if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info)
            fMagneticField->PrintMetadata();
    }
 
    RESTDebug << "Axion-field : " << fAxionField << RESTendl;
    RESTDebug << "Buffer gas : " << fBufferGas << RESTendl;
    RESTDebug << "Magnetic field : " << fMagneticField << RESTendl;
}
 
TRestEvent* TRestAxionFieldPropagationProcess::ProcessEvent(TRestEvent* evInput) {
    fAxionEvent = (TRestAxionEvent*)evInput;
 
    RESTDebug << "TRestAxionFieldPropagationProcess::ProcessEvent : " << fAxionEvent->GetID() << RESTendl;
 
    fMagneticField->SetTrack(fAxionEvent->GetPosition(), fAxionEvent->GetDirection());
 
    Double_t Ea = fAxionEvent->GetEnergy();
    Double_t ma = fAxionEvent->GetMass() * units("eV");
 
    std::pair<Double_t, Double_t> prob =
        fAxionField->GammaTransmissionFieldMapProbability(Ea, ma, fAccuracy, fNumIntervals, fQawoLevels);
 
    SetObservableValue("probability", prob.first);
    SetObservableValue("error", prob.second);
 
    Double_t transmission = 1;
    if (fBufferGas && fBufferGasAdditionalLength > 0) {
        Double_t Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea);  // cm-1
        Double_t GammaL = Gamma * fBufferGasAdditionalLength * units("cm");
        transmission = exp(-GammaL);
    }
 
    SetObservableValue("transmission", transmission);
 
    RESTDebug << " --- Process observables: " << RESTendl;
    RESTDebug << "Probability: " << prob.first << " T" << RESTendl;
    RESTDebug << "Error: " << prob.second << " T" << RESTendl;
    RESTDebug << "Transmission: " << transmission << " mm" << RESTendl;
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fAxionEvent->PrintEvent();
 
 
    return fAxionEvent;
}