TRestAxionGeneratorProcess.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 "TRestAxionGeneratorProcess.h"
using namespace std;
 
ClassImp(TRestAxionGeneratorProcess);
 
TRestAxionGeneratorProcess::TRestAxionGeneratorProcess() { Initialize(); }
 
TRestAxionGeneratorProcess::~TRestAxionGeneratorProcess() { delete fOutputAxionEvent; }
 
void TRestAxionGeneratorProcess::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fOutputAxionEvent = new TRestAxionEvent();
 
    fIsExternal = true;
    fSingleThreadOnly = true;
}
 
void TRestAxionGeneratorProcess::InitProcess() {
    RESTDebug << "Entering ... TRestAxionGeneratorProcess::InitProcess" << RESTendl;
 
    fAxionFlux = GetMetadata<TRestAxionSolarFlux>();
 
    if (fGeneratorType == "solarFlux" && fAxionFlux == nullptr) {
        if (!this->GetError()) this->SetError("The solar flux definition was not found.");
    }
 
    if (fAxionFlux) {
        fAxionFlux->Initialize();
 
        fTotalFlux = fAxionFlux->IntegrateFluxInRange(fEnergyRange);
    }
 
    if (!fRandom) {
        delete fRandom;
        fRandom = nullptr;
    }
 
    fRandom = new TRandom3(fSeed);
    fSeed = fRandom->TRandom::GetSeed();
}
 
TRestEvent* TRestAxionGeneratorProcess::ProcessEvent(TRestEvent* evInput) {
    RESTDebug << "TRestAxionGeneratorProcess::ProcessEvent : " << fCounter << RESTendl;
    fOutputAxionEvent->SetID(fCounter);
    fCounter++;
 
    TVector3 axionPosition = TVector3(0, 0, -REST_Physics::AU);
    TVector3 axionDirection = TVector3(0, 0, 1);
    Double_t energy = 1;
 
    // Random unit position
    // We avoid the use of expensive trigonometric functions
    Double_t x, y, r;
    do {
        x = 2 * (fRandom->Rndm() - 0.5);
        y = 2 * (fRandom->Rndm() - 0.5);
        r = x * x + y * y;
    } while (r > 1 || r == 0);
 
    r = TMath::Sqrt(r);
 
    std::pair<Double_t, Double_t> p;
    if (fGeneratorType == "solarFlux") {
        p = fAxionFlux->GetRandomEnergyAndRadius(fEnergyRange);
        energy = p.first;
        Double_t radius = p.second;
 
        axionPosition = TVector3(REST_Physics::solarRadius * radius * x / r,
                                 REST_Physics::solarRadius * radius * y / r, -REST_Physics::AU);
 
        axionDirection = -axionPosition.Unit();
    }
 
    if (fGeneratorType == "flat" || fGeneratorType == "plain") {
        if (fEnergyRange.X() > 0 && fEnergyRange.Y() > 0)
            energy = fRandom->Rndm() * (fEnergyRange.Y() - fEnergyRange.X()) + fEnergyRange.X();
        else if (fEnergyRange.X() > 0)
            energy = (1. - fEnergyRange.X()) * fRandom->Rndm() + fEnergyRange.X();
        else {
            RESTWarning << "Not a valid energy range was defined!" << RESTendl;
        }
    }
 
    do {
        x = 2 * (fRandom->Rndm() - 0.5);
        y = 2 * (fRandom->Rndm() - 0.5);
        r = x * x + y * y;
    } while (r > 1 || r == 0);
 
    r = TMath::Sqrt(r);
 
    axionPosition = axionPosition + TVector3(fTargetRadius * x, fTargetRadius * y, 0) + fTargetPosition;
 
    Double_t mass = fRandom->Uniform(fAxionMassRange.X(), fAxionMassRange.Y());
 
    fOutputAxionEvent->SetEnergy(energy);
    fOutputAxionEvent->SetPosition(axionPosition);
    fOutputAxionEvent->SetDirection(axionDirection);
    fOutputAxionEvent->SetMass(mass);
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug)
        fOutputAxionEvent->PrintEvent();
 
    return fOutputAxionEvent;
}
 
void TRestAxionGeneratorProcess::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    RESTMetadata << "Generator type: " << fGeneratorType << RESTendl;
    RESTMetadata << "Axion mass range: (" << fAxionMassRange.X() * units("eV") << ", "
                 << fAxionMassRange.Y() * units("eV") << ") eV" << RESTendl;
    RESTMetadata << "Target radius: " << fTargetRadius * units("cm") << " cm" << RESTendl;
    RESTMetadata << "Random seed: " << (UInt_t)fSeed << RESTendl;
    RESTMetadata << "Energy range: (" << fEnergyRange.X() << ", " << fEnergyRange.Y() << ") keV" << RESTendl;
    RESTMetadata << "Total flux: " << fTotalFlux << " cm-2 s-1" << RESTendl;
 
    RESTMetadata << "+++++++++++++++++++++++++++++++++++++++++++++++++" << RESTendl;
}