TRestAxionXrayWindow.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 "TRestAxionXrayWindow.h"
using namespace std;
 
#include "TRestSystemOfUnits.h"
using namespace REST_Units;
 
ClassImp(TRestAxionXrayWindow);
 
TRestAxionXrayWindow::TRestAxionXrayWindow() : TRestMetadata() { Initialize(); }
 
TRestAxionXrayWindow::TRestAxionXrayWindow(const char* cfgFileName, string name)
    : TRestMetadata(cfgFileName) {
    RESTDebug << "Entering TRestAxionXrayWindow constructor( cfgFileName, name )" << RESTendl;
 
    Initialize();
 
    LoadConfigFromFile(fConfigFileName, name);
}
 
TRestAxionXrayWindow::~TRestAxionXrayWindow() {}
 
void TRestAxionXrayWindow::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fEnergy.clear();
    fTransmission.clear();
 
    /*
if (fWindowType != "foil" && fPatternGap == 0) {
    RESTError << "TRestAxionXrayWindow::Initialize. fPatternGap cannot be zero!" << RESTendl;
    fPatternGap = 1;
}
 
if (fPatternGap < 0) fPatternGap = -fPatternGap;
    */
}
 
void TRestAxionXrayWindow::ReadMaterial() {
    std::string materialFileName = SearchFile(fMaterial + ".sol");
 
    RESTDebug << "TRestAxionXrayWindow::ReadMaterial. Reading file : " << materialFileName << RESTendl;
 
    if (!TRestTools::fileExists(materialFileName)) {
        RESTError << "TRestAxionXrayWindow::ReadMaterial( )" << RESTendl;
        RESTError << "Material file not found : " << materialFileName << RESTendl;
        exit(1);
    }
 
    FILE* fin = fopen(materialFileName.c_str(), "rt");
 
    fEnergy.clear();
    fTransmission.clear();
    double en, value;
    while (fscanf(fin, "%lf\t%lf\n", &en, &value) != EOF) {
        RESTDebug << "Energy : " << en << "eV -- Abs : " << value << RESTendl;
 
        fEnergy.push_back(en / 1000.);
        fTransmission.push_back(TMath::Power(value, fThickness * units("um")));
    }
 
    RESTDebug << "Items read : " << fEnergy.size() << RESTendl;
 
    fclose(fin);
}
 
Double_t TRestAxionXrayWindow::GetTransmission(Double_t energy, Double_t x, Double_t y) {
    if (fMaterial != "Vacuum" && fEnergy.size() == 0) ReadMaterial();
 
    if (fMask && x * x + y * y > fMask->GetMaskRadius() * fMask->GetMaskRadius()) return 0;
 
    if (fMask && !fMask->HitsPattern(x, y)) return 1.;
 
    // This line must be after all previous cuts to make sure mask limits do an effect
    if (fMaterial == "Vacuum") return 1;
 
    Double_t energyIndex = GetEnergyIndex(energy);
 
    if (energyIndex < 0) {
        RESTWarning << "Energy : " << energy << " keV is out of range!" << RESTendl;
        return 0;
    }
 
    // Transmission
    double y2 = fTransmission[energyIndex + 1];
    double y1 = fTransmission[energyIndex];
 
    // Normalized field
    double x2 = fEnergy[energyIndex + 1];
    double x1 = fEnergy[energyIndex];
 
    double m = (y2 - y1) / (x2 - x1);
    double n = y1 - m * x1;
 
    if (m * energy + n < 0) {
        RESTError << "TRestAxionXrayWindow::GetAbsorptionCoefficient. Negative coefficient!" << RESTendl;
        cout << "y2 : " << y2 << " y1 : " << y1 << endl;
        cout << "x2 : " << x2 << " x1 : " << x1 << endl;
        cout << "m : " << m << " n : " << n << endl;
        cout << "E : " << energy << " bin : " << energyIndex << endl;
        return 0;
    }
 
    return (m * energy + n);
}
 
Bool_t TRestAxionXrayWindow::HitsPattern(Double_t x, Double_t y) {
    if (fMask) return fMask->HitsPattern(x, y);
 
    /*
if (fWindowType == "stripped") {
    Double_t xEval = fPatternWidth / 2. + x - fPatternOffset;
 
    if (xEval > 0) {
        while (xEval > fPatternGap) xEval -= fPatternGap;
    } else {
        while (xEval < 0) xEval += fPatternGap;
    }
 
    if (xEval > fPatternWidth) {
        return false;
    }
} else if (fWindowType == "grid") {
    Double_t xEval = fPatternWidth / 2. + x - fPatternOffset;
 
    if (xEval > 0) {
        while (xEval > fPatternGap) xEval -= fPatternGap;
    } else {
        while (xEval < 0) xEval += fPatternGap;
    }
 
    if (xEval < fPatternWidth) return true;
 
    Double_t yEval = fPatternWidth / 2. + y - fPatternOffset;
 
    if (yEval > 0) {
        while (yEval > fPatternGap) yEval -= fPatternGap;
    } else {
        while (yEval < 0) yEval += fPatternGap;
    }
 
    if (yEval < fPatternWidth) return true;
    return false;
}
    */
 
    return true;
}
 
Int_t TRestAxionXrayWindow::GetEnergyIndex(Double_t energy) {
    for (unsigned int n = 0; n < fEnergy.size(); n++)
        if (energy < fEnergy[n]) return n - 1;
    return -1;
}
 
void TRestAxionXrayWindow::PrintTransmissionData() {
    if (fEnergy.size() == 0) ReadMaterial();
    for (unsigned int n = 0; n < fEnergy.size(); n++)
        cout << "Energy : " << fEnergy[n] << " Transmission : " << fTransmission[n] << endl;
}
 
void TRestAxionXrayWindow::InitFromConfigFile() {
    TRestMetadata::InitFromConfigFile();
 
    if (fMask) {
        delete fMask;
        fMask = nullptr;
    }
    fMask = (TRestPatternMask*)this->InstantiateChildMetadata(0, "Mask");
 
    if (!fMask) {
        RESTWarning << "TRestAxionXrayWindow. Name : " << this->GetName() << RESTendl;
        RESTWarning << "No mask pattern was defined for the X-ray window!" << RESTendl;
    }
}
 
void TRestAxionXrayWindow::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    RESTMetadata << "Thickness: " << fThickness * units("um") << " um" << RESTendl;
    RESTMetadata << "Material: " << fMaterial << RESTendl;
    RESTMetadata << "----" << RESTendl;
    if (fMask) {
        fMask->Print();
    } else {
        RESTMetadata << " - Pattern type: foil" << RESTendl;
    }
 
    RESTMetadata << "+++++++++++++++++++++++++++++++++++++++++++++++++" << RESTendl;
}