TRestAxionSolarQCDFlux.cxx

/******************** REST disclaimer ***********************************
 * 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 "TRestAxionSolarQCDFlux.h"
using namespace std;
 
ClassImp(TRestAxionSolarQCDFlux);
 
TRestAxionSolarQCDFlux::TRestAxionSolarQCDFlux() : TRestAxionSolarFlux() {}
 
TRestAxionSolarQCDFlux::TRestAxionSolarQCDFlux(const char* cfgFileName, string name)
    : TRestAxionSolarFlux(cfgFileName) {
    LoadConfigFromFile(fConfigFileName, name);
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata();
}
 
TRestAxionSolarQCDFlux::~TRestAxionSolarQCDFlux() {}
 
Bool_t TRestAxionSolarQCDFlux::LoadTables() {
    if (fFluxDataFile == "" && fFluxSptFile == "") return false;
 
    if (TRestTools::GetFileNameExtension(fFluxDataFile) == "flux") {
        ReadFluxFile();
    } else {
        LoadContinuumFluxTable();
        LoadMonoChromaticFluxTable();
    }
 
    IntegrateSolarFluxes();
 
    return true;
}
 
void TRestAxionSolarQCDFlux::LoadContinuumFluxTable() {
    if (fFluxDataFile == "") {
        RESTDebug << "TRestAxionSolarflux::LoadContinuumFluxTable. No solar flux table was defined"
                  << RESTendl;
        return;
    }
 
    string fullPathName = SearchFile((string)fFluxDataFile);
 
    RESTDebug << "Loading table from file : " << RESTendl;
    RESTDebug << "File : " << fullPathName << RESTendl;
 
    std::vector<std::vector<Float_t>> fluxTable;
    if (TRestTools::GetFileNameExtension(fFluxDataFile) == "dat") {
        std::vector<std::vector<Double_t>> doubleTable;
        if (!TRestTools::ReadASCIITable(fullPathName, doubleTable)) {
            RESTError << "TRestAxionSolarQCDFlux::LoadContinuumFluxTable. " << RESTendl;
            RESTError << "Could not read solar flux table : " << fFluxDataFile << RESTendl;
            return;
        }
        for (const auto& row : doubleTable) {
            std::vector<Float_t> floatVec(row.begin(), row.end());
            fluxTable.push_back(floatVec);
        }
    } else if (TRestTools::IsBinaryFile(fFluxDataFile))
        TRestTools::ReadBinaryTable(fullPathName, fluxTable);
    else {
        fluxTable.clear();
        RESTError << "Filename extension was not recognized!" << RESTendl;
        RESTError << "Solar flux table will not be populated" << RESTendl;
        RESTError << "Filename extension: " << TRestTools::GetFileNameExtension(fFluxDataFile) << RESTendl;
    }
 
    if (fluxTable.size() != 100 && fluxTable[0].size() != 200) {
        fluxTable.clear();
        RESTError << "LoadContinuumFluxTable. The table does not contain the right number of rows or columns"
                  << RESTendl;
        RESTError << "Table will not be populated" << RESTendl;
    }
 
    for (unsigned int n = 0; n < fluxTable.size(); n++) {
        TH1F* h = new TH1F(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), "", 200, 0, 20);
        for (unsigned int m = 0; m < fluxTable[n].size(); m++) h->SetBinContent(m + 1, fluxTable[n][m]);
        fFluxTable.push_back(h);
    }
}
 
void TRestAxionSolarQCDFlux::LoadMonoChromaticFluxTable() {
    if (fFluxSptFile == "") {
        RESTDebug << "TRestAxionSolarflux::LoadMonoChromaticFluxTable. No solar flux monochromatic table was "
                     "defined"
                  << RESTendl;
        return;
    }
 
    string fullPathName = SearchFile((string)fFluxSptFile);
 
    RESTDebug << "Loading monochromatic lines from file : " << RESTendl;
    RESTDebug << "File : " << fullPathName << RESTendl;
 
    std::vector<std::vector<Double_t>> doubleTable;
    if (!TRestTools::ReadASCIITable(fullPathName, doubleTable)) {
        RESTError << "TRestAxionSolarQCDFlux::LoadMonoChromaticFluxTable." << RESTendl;
        RESTError << "Could not read solar flux table : " << fFluxSptFile << RESTendl;
        return;
    }
 
    std::vector<std::vector<Float_t>> asciiTable;
    for (const auto& row : doubleTable) {
        std::vector<Float_t> floatVec(row.begin(), row.end());
        asciiTable.push_back(floatVec);
    }
 
    fFluxLines.clear();
 
    if (asciiTable.size() != 101) {
        RESTError << "LoadMonoChromaticFluxTable. The table does not contain the right number of rows"
                  << RESTendl;
        RESTError << "Table will not be populated" << RESTendl;
        return;
    }
 
    for (unsigned int en = 0; en < asciiTable[0].size(); en++) {
        Float_t energy = asciiTable[0][en];
        TH1F* h = new TH1F(Form("%s_MonochromeFluxAtEnergy%6.4lf", GetName(), energy), "", 100, 0, 1);
        for (unsigned int r = 1; r < asciiTable.size(); r++) h->SetBinContent(r, asciiTable[r][en]);
        fFluxLines[energy] = h;
    }
}
 
void TRestAxionSolarQCDFlux::ReadFluxFile() {
    if (fBinSize <= 0) {
        RESTError << "TRestAxionSolarflux::ReadFluxFile. Energy bin size of .flux file must be specified."
                  << RESTendl;
        RESTError << "Please, define binSize parameter in eV." << RESTendl;
        return;
    }
 
    if (fPeakSigma <= 0) {
        RESTWarning << "TRestAxionSolarflux::ReadFluxFile. Peak sigma must be specified to generate "
                       "monochromatic spectrum."
                    << RESTendl;
        RESTWarning << "Only continuum table will be generated. If this was intentional, please, ignore this "
                       "RESTWarning."
                    << RESTendl;
        return;
    }
 
    string fullPathName = SearchFile((string)fFluxDataFile);
 
    RESTDebug << "Loading flux table ...  " << RESTendl;
    RESTDebug << "File : " << fullPathName << RESTendl;
    std::vector<std::vector<Double_t>> fluxData;
    TRestTools::ReadASCIITable(fullPathName, fluxData, 3);
 
    RESTDebug << "Table loaded" << RESTendl;
    TH2F* originalHist = new TH2F("FullTable", "", 100, 0., 1., (Int_t)(20. / fBinSize), 0., 20.);
    TH2F* continuumHist = new TH2F("ContinuumTable", "", 100, 0., 1., (Int_t)(20. / fBinSize), 0., 20.);
    TH2F* spectrumHist = new TH2F("LinesTable", "", 100, 0., 1., (Int_t)(20. / fBinSize), 0., 20.);
 
    Double_t fluxBinSize = TRestTools::GetLowestIncreaseFromTable(fluxData, 1);
 
    for (const auto& data : fluxData) {
        Float_t r = 0.005 + data[0];
        Float_t en = data[1] - 0.005;
        Float_t flux = data[2] * fluxBinSize;  // flux in cm-2 s-1 bin-1
 
        originalHist->Fill(r, en, (Float_t)flux);
        continuumHist->Fill(r, en, (Float_t)flux);
    }
    RESTDebug << "Histograms filled" << RESTendl;
 
    Int_t peaks = 0;
    do {
        peaks = 0;
        // We just identify pronounced peaks, not smoothed gaussians!
        const int smearPoints = (Int_t)(5 / (fBinSize * 100));
        const int excludePoints = smearPoints / 5;
        for (const auto& data : fluxData) {
            Float_t r = 0.005 + data[0];
            Float_t en = data[1] - 0.005;
 
            Int_t binR = continuumHist->GetXaxis()->FindBin(r);
            Int_t binE = continuumHist->GetYaxis()->FindBin(en);
 
            Double_t avgFlux = 0;
            Int_t n = 0;
            for (int e = binE - smearPoints; e <= binE + smearPoints; e++) {
                if (e < 1 || (e < binE + excludePoints && e > binE - excludePoints)) continue;
                n++;
                avgFlux += continuumHist->GetBinContent(binR, e);
            }
            avgFlux /= n;
 
            Float_t targetBinFlux = continuumHist->GetBinContent(binR, binE);
            Float_t thrFlux = avgFlux + fPeakSigma * TMath::Sqrt(avgFlux);
            if (targetBinFlux > 0 && targetBinFlux > thrFlux) {
                continuumHist->SetBinContent(binR, binE, avgFlux);
                peaks++;
            }
        }
    } while (peaks > 0);
 
    for (int n = 0; n < originalHist->GetNbinsX(); n++)
        for (int m = 0; m < originalHist->GetNbinsY(); m++) {
            Float_t orig = originalHist->GetBinContent(n + 1, m + 1);
            Float_t cont = continuumHist->GetBinContent(n + 1, m + 1);
 
            spectrumHist->SetBinContent(n + 1, m + 1, orig - cont);
        }
 
    continuumHist->Rebin2D(1, (Int_t)(0.1 / fBinSize));  // cm-2 s-1 (100eV)-1
    continuumHist->Scale(10);                            // cm-2 s-1 keV-1
    // It could be over here if we would use directly a TH2D
 
    fFluxTable.clear();
    for (int n = 0; n < continuumHist->GetNbinsX(); n++) {
        TH1F* hc =
            (TH1F*)continuumHist->ProjectionY(Form("%s_ContinuumFluxAtRadius%d", GetName(), n), n + 1, n + 1);
        fFluxTable.push_back(hc);
    }
 
    fFluxLines.clear();
    for (int n = 0; n < spectrumHist->GetNbinsY(); n++) {
        if (spectrumHist->ProjectionX("", n + 1, n + 1)->Integral() > 0) {
            Double_t energy = spectrumHist->ProjectionY()->GetBinCenter(n + 1);
            TH1F* hm = (TH1F*)spectrumHist->ProjectionX(
                Form("%s_MonochromeFluxAtEnergy%5.3lf", GetName(), energy), n + 1, n + 1);
            fFluxLines[energy] = hm;
        }
    }
 
    cout << "Number of peaks identified: " << fFluxLines.size() << endl;
}
 
TH1F* TRestAxionSolarQCDFlux::GetContinuumSpectrum() {
    if (fContinuumHist != nullptr) {
        delete fContinuumHist;
        fContinuumHist = nullptr;
    }
 
    fContinuumHist = new TH1F("ContinuumHist", "", 200, 0, 20);
    for (const auto& x : fFluxTable) {
        fContinuumHist->Add(x);
    }
 
    fContinuumHist->SetStats(0);
    fContinuumHist->GetXaxis()->SetTitle("Energy [keV]");
    fContinuumHist->GetXaxis()->SetTitleSize(0.05);
    fContinuumHist->GetXaxis()->SetLabelSize(0.05);
    fContinuumHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 keV-1]");
    fContinuumHist->GetYaxis()->SetTitleSize(0.05);
    fContinuumHist->GetYaxis()->SetLabelSize(0.05);
 
    return fContinuumHist;
}
 
TH1F* TRestAxionSolarQCDFlux::GetMonochromaticSpectrum() {
    if (fMonoHist != nullptr) {
        delete fMonoHist;
        fMonoHist = nullptr;
    }
 
    fMonoHist = new TH1F("MonochromaticHist", "", 20000, 0, 20);
    for (const auto& x : fFluxLines) {
        fMonoHist->Fill(x.first, x.second->Integral());  // cm-2 s-1 eV-1
    }
 
    fMonoHist->SetStats(0);
    fMonoHist->GetXaxis()->SetTitle("Energy [keV]");
    fMonoHist->GetXaxis()->SetTitleSize(0.05);
    fMonoHist->GetXaxis()->SetLabelSize(0.05);
    fMonoHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 eV-1]");
    fMonoHist->GetYaxis()->SetTitleSize(0.05);
    fMonoHist->GetYaxis()->SetLabelSize(0.05);
 
    return fMonoHist;
}
 
TH1F* TRestAxionSolarQCDFlux::GetTotalSpectrum() {
    TH1F* hm = GetMonochromaticSpectrum();
    TH1F* hc = GetContinuumSpectrum();
 
    if (fTotalHist != nullptr) {
        delete fTotalHist;
        fTotalHist = nullptr;
    }
 
    fTotalHist = new TH1F("fTotalHist", "", 20000, 0, 20);
    for (int n = 0; n < hc->GetNbinsX(); n++) {
        for (int m = 0; m < 100; m++) {
            fTotalHist->SetBinContent(n * 100 + 1 + m, hc->GetBinContent(n + 1));
        }
    }
 
    for (int n = 0; n < hm->GetNbinsX(); n++)
        // 1e-2 is the renormalization from 20000 bins to 200 bins
        fTotalHist->SetBinContent(n + 1, fTotalHist->GetBinContent(n + 1) + 100 * hm->GetBinContent(n + 1));
 
    fTotalHist->SetStats(0);
    fTotalHist->GetXaxis()->SetTitle("Energy [keV]");
    fTotalHist->GetXaxis()->SetTitleSize(0.05);
    fTotalHist->GetXaxis()->SetLabelSize(0.05);
    fTotalHist->GetYaxis()->SetTitle("Flux [cm-2 s-1 keV-1]");
    fTotalHist->GetYaxis()->SetTitleSize(0.05);
    fTotalHist->GetYaxis()->SetLabelSize(0.05);
 
    return fTotalHist;
}
 
TCanvas* TRestAxionSolarQCDFlux::DrawSolarFlux() {
    if (fCanvas != nullptr) {
        delete fCanvas;
        fCanvas = nullptr;
    }
    fCanvas = new TCanvas("canv", "This is the canvas title", 1200, 500);
    fCanvas->Draw();
 
    TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97);
    pad1->Divide(2, 1);
    pad1->Draw();
 
    pad1->cd(1);
    pad1->cd(1)->SetLogy();
    pad1->cd(1)->SetRightMargin(0.09);
    pad1->cd(1)->SetLeftMargin(0.15);
    pad1->cd(1)->SetBottomMargin(0.15);
 
    TH1F* ht = GetTotalSpectrum();
    ht->SetLineColor(kBlack);
    ht->SetFillStyle(4050);
    ht->SetFillColor(kBlue - 10);
 
    TH1F* hm = GetMonochromaticSpectrum();
    hm->SetLineColor(kBlack);
    hm->Scale(100);  // renormalizing per 100eV-1
 
    ht->Draw("hist");
    hm->Draw("hist same");
 
    pad1->cd(2);
    pad1->cd(2)->SetRightMargin(0.09);
    pad1->cd(2)->SetLeftMargin(0.15);
    pad1->cd(2)->SetBottomMargin(0.15);
 
    ht->Draw("hist");
    hm->Draw("hist same");
 
    return fCanvas;
}
 
void TRestAxionSolarQCDFlux::IntegrateSolarFluxes() {
    fFluxLineIntegrals.clear();
    fTotalMonochromaticFlux = 0;
 
    for (const auto& line : fFluxLines) {
        fTotalMonochromaticFlux += line.second->Integral();
        fFluxLineIntegrals.push_back(fTotalMonochromaticFlux);
    }
 
    for (unsigned int n = 0; n < fFluxLineIntegrals.size(); n++)
        fFluxLineIntegrals[n] /= fTotalMonochromaticFlux;
 
    fTotalContinuumFlux = 0.0;
    for (unsigned int n = 0; n < fFluxTable.size(); n++) {
        fTotalContinuumFlux += fFluxTable[n]->Integral() * 0.1;  // We integrate in 100eV steps
        fFluxTableIntegrals.push_back(fTotalContinuumFlux);
    }
 
    for (unsigned int n = 0; n < fFluxTableIntegrals.size(); n++)
        fFluxTableIntegrals[n] /= fTotalContinuumFlux;
 
    fFluxRatio = fTotalMonochromaticFlux / (fTotalContinuumFlux + fTotalMonochromaticFlux);
}
 
Double_t TRestAxionSolarQCDFlux::IntegrateFluxInRange(TVector2 eRange, Double_t mass) {
    if (eRange.X() == -1 && eRange.Y() == -1) {
        if (GetTotalFlux() == 0) IntegrateSolarFluxes();
        return GetTotalFlux();
    }
 
    Double_t flux = 0;
    for (const auto& line : fFluxLines)
        if (line.first > eRange.X() && line.first < eRange.Y()) flux += line.second->Integral();
 
    fTotalContinuumFlux = 0.0;
    for (unsigned int n = 0; n < fFluxTable.size(); n++) {
        flux += fFluxTable[n]->Integral(fFluxTable[n]->FindFixBin(eRange.X()),
                                        fFluxTable[n]->FindFixBin(eRange.Y())) *
                0.1;  // We integrate in 100eV steps
    }
 
    return flux;
}
 
std::pair<Double_t, Double_t> TRestAxionSolarQCDFlux::GetRandomEnergyAndRadius(TVector2 eRange,
                                                                               Double_t mass) {
    std::pair<Double_t, Double_t> result = {0, 0};
    if (!AreTablesLoaded()) return result;
    Double_t rnd = fRandom->Rndm();
    if (fTotalMonochromaticFlux == 0 || fRandom->Rndm() > fFluxRatio) {
        // Continuum
        for (unsigned int r = 0; r < fFluxTableIntegrals.size(); r++) {
            if (rnd < fFluxTableIntegrals[r]) {
                Double_t energy = fFluxTable[r]->GetRandom();
                if (eRange.X() != -1 && eRange.Y() != -1) {
                    if (energy < eRange.X() || energy > eRange.Y()) return GetRandomEnergyAndRadius(eRange);
                }
                Double_t radius = ((Double_t)r + fRandom->Rndm()) * 0.01;
                std::pair<Double_t, Double_t> p = {energy, radius};
                return p;
            }
        }
    } else {
        // Monochromatic
        int n = 0;
        for (const auto& line : fFluxLines) {
            if (rnd < fFluxLineIntegrals[n]) {
                std::pair<Double_t, Double_t> p = {line.first, line.second->GetRandom()};
                return p;
            }
            n++;
        }
    }
    return result;
}
 
void TRestAxionSolarQCDFlux::PrintContinuumSolarTable() {
    cout << "Continuum solar flux table: " << endl;
    cout << "--------------------------- " << endl;
    for (unsigned int n = 0; n < fFluxTable.size(); n++) {
        for (int m = 0; m < fFluxTable[n]->GetNbinsX(); m++)
            cout << fFluxTable[n]->GetBinContent(m + 1) << "\t";
        cout << endl;
        cout << endl;
    }
    cout << endl;
}
 
void TRestAxionSolarQCDFlux::PrintIntegratedRingFlux() {
    cout << "Integrated solar flux per solar ring: " << endl;
    cout << "--------------------------- " << endl;
    /*
    for (int n = 0; n < fFluxPerRadius.size(); n++)
    cout << "n : " << n << " flux : " << fFluxPerRadius[n] << endl;
    cout << endl;
    */
}
 
void TRestAxionSolarQCDFlux::PrintMonoChromaticFlux() {
    //   cout << "Number of monochromatic lines: " << fFluxPerRadius.size() << endl;
    cout << "+++++++++++++++++++++++++++++++++++" << endl;
    for (auto const& line : fFluxLines) {
        cout << "Energy : " << line.first << " keV" << endl;
        cout << "-----------------" << endl;
        for (int n = 0; n < line.second->GetNbinsX(); n++)
            cout << "R : " << line.second->GetBinCenter(n + 1)
                 << " flux : " << line.second->GetBinContent(n + 1) << " cm-2 s-1" << endl;
    }
}
 
void TRestAxionSolarQCDFlux::PrintMetadata() {
    TRestAxionSolarFlux::PrintMetadata();
 
    if (fFluxDataFile != "")
        RESTMetadata << " - Solar axion flux datafile (continuum) : " << fFluxDataFile << RESTendl;
    if (fFluxSptFile != "")
        RESTMetadata << " - Solar axion flux datafile (monochromatic) : " << fFluxSptFile << RESTendl;
    RESTMetadata << "-------" << RESTendl;
    RESTMetadata << " - Total monochromatic flux : " << fTotalMonochromaticFlux << " cm-2 s-1" << RESTendl;
    RESTMetadata << " - Total continuum flux : " << fTotalContinuumFlux << " cm-2 s-1" << RESTendl;
    RESTMetadata << "++++++++++++++++++" << RESTendl;
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) {
        PrintContinuumSolarTable();
        PrintMonoChromaticFlux();
        PrintIntegratedRingFlux();
    }
}
 
void TRestAxionSolarQCDFlux::ExportTables(Bool_t ascii) {
    string rootFilename = TRestTools::GetFileNameRoot(fFluxDataFile);
 
    string path = REST_USER_PATH + "/export/";
 
    if (!TRestTools::fileExists(path)) {
        std::cout << "Creating path: " << path << std::endl;
        int z = system(("mkdir -p " + path).c_str());
        if (z != 0) RESTError << "Could not create directory " << path << RESTendl;
    }
 
    if (fFluxTable.size() > 0) {
        std::vector<std::vector<Float_t>> table;
        for (const auto& x : fFluxTable) {
            std::vector<Float_t> row;
            for (int n = 0; n < x->GetNbinsX(); n++) row.push_back(x->GetBinContent(n + 1));
 
            table.push_back(row);
        }
 
        if (!ascii)
            TRestTools::ExportBinaryTable(path + "/" + rootFilename + ".N200f", table);
        else
            TRestTools::ExportASCIITable(path + "/" + rootFilename + ".dat", table);
    }
 
    if (fFluxLines.size() > 0) {
        std::vector<std::vector<Float_t>> table;
        for (const auto& x : fFluxLines) {
            std::vector<Float_t> row;
            row.push_back(x.first);
            for (int n = 0; n < x.second->GetNbinsX(); n++) row.push_back(x.second->GetBinContent(n + 1));
 
            table.push_back(row);
        }
 
        TRestTools::TransposeTable(table);
 
        TRestTools::ExportASCIITable(path + "/" + rootFilename + ".spt", table);
    }
}