src/TRestWimpSensitivity.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 "TRestWimpSensitivity.h"
#include "TRestWimpUtils.h"
 
ClassImp(TRestWimpSensitivity);
 
using namespace std;
 
TRestWimpSensitivity::TRestWimpSensitivity(const char* configFilename, const string& name)
    : TRestMetadata(configFilename) {
    Initialize();
 
    LoadConfigFromFile(fConfigFileName, name);
}
TRestWimpSensitivity::~TRestWimpSensitivity() {}
 
void TRestWimpSensitivity::Initialize() {
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
}
 
void TRestWimpSensitivity::InitFromConfigFile() {
    this->Initialize();
    TRestMetadata::InitFromConfigFile();
    ReadNuclei();
}
 
void TRestWimpSensitivity::ReadNuclei() {
    fNuclei.clear();
 
    TiXmlElement* ElementDef = GetElement("addElement");
    while (ElementDef) {
        TRestWimpNucleus nucleus;
        nucleus.fNucleusName = GetFieldValue("nucleusName", ElementDef);
        nucleus.fAnum = StringToDouble(GetFieldValue("anum", ElementDef));
        nucleus.fZnum = StringToInteger(GetFieldValue("znum", ElementDef));
        nucleus.fAbundance = StringToDouble(GetFieldValue("abundance", ElementDef));
        fNuclei.emplace_back(nucleus);
        ElementDef = GetNextElement(ElementDef);
    }
}
 
std::map<std::string, TH1D*> TRestWimpSensitivity::GetRecoilSpectra(const double wimpMass,
                                                                    const double crossSection) {
    std::map<std::string, TH1D*> recoilMap;
 
    const int nBins = (fEnergySpectra.Y() - fEnergySpectra.X()) / fEnergySpectraStep;
 
    for (auto& nucl : fNuclei) {
        std::string histName = "RecoilSpc_" + std::string(nucl.fNucleusName);
        TH1D* recoilSpc =
            new TH1D(histName.c_str(), histName.c_str(), nBins, fEnergySpectra.X(), fEnergySpectra.Y());
        for (int i = 1; i < recoilSpc->GetNbinsX(); i++) {
            const double recoilEnergy = recoilSpc->GetBinCenter(i);
            const double recoilRate =
                TRestWimpUtils::GetRecoilRate(wimpMass, crossSection, recoilEnergy, nucl.fAnum, fLabVelocity,
                                              fRmsVelocity, fEscapeVelocity, fWimpDensity, nucl.fAbundance);
            recoilSpc->SetBinContent(i, recoilRate);
        }
        recoilMap[std::string(nucl.fNucleusName)] = recoilSpc;
    }
 
    return recoilMap;
}
 
const Double_t TRestWimpSensitivity::GetSensitivity(const double wimpMass) {
    if (fNuclei.empty()) {
        RESTError << "Please add at least one element to the rml file" << RESTendl;
    }
 
    if (fUseQuenchingFactor) CalculateQuenchingFactor();
 
    const int nBins = (fEnergySpectra.Y() - fEnergySpectra.X()) / fEnergySpectraStep;
 
    TH1D recoilSpc("recoilSpc", "recoilSpc", nBins, fEnergySpectra.X(), fEnergySpectra.Y());
 
    double bckCounts = 0;
 
    for (int i = 1; i < recoilSpc.GetNbinsX(); i++) {
        const double en = recoilSpc.GetBinCenter(i);
        if (en < fEnergyRange.X() || en > fEnergyRange.Y()) continue;
        bckCounts += fBackground * fEnergySpectraStep;
    }
    bckCounts *= fExposure;
 
    double signalCounts = 0, prob = 0;
 
    do {
        prob = 0;
        for (int mu = signalCounts; mu < (signalCounts + bckCounts + 10000); mu++) {
            if (bckCounts <= 1.e3)
                prob += TMath::Poisson(mu + bckCounts, bckCounts);
            else if (bckCounts > 1.e3)
                prob += TMath::Gaus(mu + bckCounts, bckCounts, TMath::Sqrt(bckCounts), true);
        }
        signalCounts++;
    } while (fabs(prob - 0.1) > 0.01 && signalCounts < 1E6);
 
    double nMeas = 0;
    const double crossSection = 1E-45;
 
    for (auto& nucl : fNuclei) {
        recoilSpc.Reset();
 
        for (int i = 1; i < recoilSpc.GetNbinsX(); i++) {
            const double recoilEnergy = recoilSpc.GetBinCenter(i);
            const double recoilRate =
                TRestWimpUtils::GetRecoilRate(wimpMass, crossSection, recoilEnergy, nucl.fAnum, fLabVelocity,
                                              fRmsVelocity, fEscapeVelocity, fWimpDensity, nucl.fAbundance);
            recoilSpc.SetBinContent(i, recoilRate);
        }
 
        for (int i = 1; i < recoilSpc.GetNbinsX(); i++) {
            double recoilEnergy = recoilSpc.GetBinCenter(i);
            // const double recoilRate = recoilSpc.GetBinContent(i);
            if (fUseQuenchingFactor)
                recoilEnergy *= quenchingFactor[std::string(nucl.fNucleusName)]->GetBinContent(i);
 
            if (recoilEnergy < fEnergyRange.X() || recoilEnergy > fEnergyRange.Y()) continue;
            nMeas += recoilSpc.GetBinContent(i) * fEnergySpectraStep;
        }
    }
 
    if (nMeas == 0) return 0;
 
    const double sensitivity = signalCounts * 1E-45 / (nMeas * fExposure);
 
    return sensitivity;
}
 
void TRestWimpSensitivity::CalculateQuenchingFactor() {
    if (!quenchingFactor.empty()) return;
 
    std::cout << "Calculating quenching factor " << std::endl;
 
    const int nBins = (fEnergySpectra.Y() - fEnergySpectra.X()) / fEnergySpectraStep;
 
    for (auto& nucl : fNuclei) {
        std::string histName = "QF_" + std::string(nucl.fNucleusName);
        TH1D* QF =
            new TH1D(histName.c_str(), histName.c_str(), nBins, fEnergySpectra.X(), fEnergySpectra.Y());
        for (int i = 1; i < QF->GetNbinsX(); i++) {
            const double recoilEnergy = QF->GetBinCenter(i);
            const double qF = TRestWimpUtils::GetQuenchingFactor(recoilEnergy, nucl.fAnum, nucl.fZnum);
            if (!isnan(qF) && qF > 0)
                QF->SetBinContent(i, 1. / qF);
            else
                QF->SetBinContent(i, 0);
        }
        quenchingFactor[std::string(nucl.fNucleusName)] = QF;
    }
}
 
const std::string TRestWimpSensitivity::BuildOutputFileName(const std::string& extension) {
    std::stringstream ss;
    ss << "WimpSensitivity_";
 
    for (auto& nucl : fNuclei) ss << nucl.fNucleusName << "_" << nucl.fAbundance << "_";
 
    ss << "WD_" << fWimpDensity << "_";
    ss << "Vel_" << fLabVelocity << "_" << fRmsVelocity << "_" << fEscapeVelocity << "_";
    ss << "Bck_" << fBackground << "_";
    ss << "RecEn_" << fEnergySpectra.X() << "_" << fEnergySpectra.Y() << "_" << fEnergySpectraStep << "_";
    ss << "EnRange_" << fEnergyRange.X() << "_" << fEnergyRange.Y() << "_";
 
    if (fUseQuenchingFactor)
        ss << "usingQF";
    else
        ss << "noQF";
 
    ss << extension;
 
    std::cout << "Output File " << ss.str() << std::endl;
 
    return ss.str();
}
 
std::map<std::string, TH1D*> TRestWimpSensitivity::GetFormFactor() {
    std::map<std::string, TH1D*> formFactor;
 
    const int nBins = (fEnergySpectra.Y() - fEnergySpectra.X()) / fEnergySpectraStep;
 
    for (auto& nucl : fNuclei) {
        std::string histName = "FF_" + std::string(nucl.fNucleusName);
        TH1D* FF =
            new TH1D(histName.c_str(), histName.c_str(), nBins, fEnergySpectra.X(), fEnergySpectra.Y());
        for (int i = 1; i < FF->GetNbinsX(); i++) {
            const double recoilEnergy = FF->GetBinCenter(i);
            const double helmFF = TRestWimpUtils::GetHelmFormFactor(recoilEnergy, nucl.fAnum);
            FF->SetBinContent(i, helmFF * helmFF);
        }
        formFactor[std::string(nucl.fNucleusName)] = FF;
    }
 
    return formFactor;
}
 
void TRestWimpSensitivity::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    for (auto& nucl : fNuclei) nucl.PrintNucleus();
 
    RESTMetadata << "WimpDensity: " << fWimpDensity << " GeV/cm3" << RESTendl;
    RESTMetadata << "WimpVelocity; VLab: " << fLabVelocity << " VRMS: " << fRmsVelocity
                 << " VEscape: " << fEscapeVelocity << " km/s" << RESTendl;
    RESTMetadata << "Exposure: " << fExposure << " kg*day" << RESTendl;
    RESTMetadata << "Background Level: " << fBackground << " c/keV/day" << RESTendl;
    RESTMetadata << "Recoil energy range: (" << fEnergySpectra.X() << ", " << fEnergySpectra.Y()
                 << ") Step: " << fEnergySpectraStep << " keV" << RESTendl;
    RESTMetadata << "Sensitivity energy range: (" << fEnergyRange.X() << ", " << fEnergyRange.Y() << ") keV"
                 << RESTendl;
    RESTMetadata << "Use quenching factor: " << fUseQuenchingFactor << RESTendl;
    RESTMetadata << "+++++" << RESTendl;
}