Few updates adding plots and minor calculations for draft paper

README.md

@@ -69,3 +69,21 @@ This repository makes use of the following published codes:
 - K. Altenmuller et al, REST-for-Physics, a ROOT-based framework for event oriented data analysis and combined Monte Carlo response, [Computer Physics Communications 273, April 2022, 108281](https://doi.org/10.1016/j.cpc.2021.108281).
 - S.Hoof, J.Jaeckel, T.J.Lennert, Quantifying uncertainties in the solar axion flux and their impact on determining axion model parameters, [JCAP09(2021)006](https://doi.org/10.1088/1475-7516/2021/09/006).
 - T.Kittelmann, E.Klinkby, E.B.Knudsen, P.Willendrup, X.X.Cai, K.Kanaki, Monte Carlo Particle Lists: MCPL, Computer Physics Communications 218 (2017) 17–42](http://dx.doi.org/10.17632/cby92vsv5g.1).
+
+### Gallery
+
+<p align="center">
+
+<img src="./images/paper/cover1.png" alt="cover1" style="width:450px;" class="center" />
+<img src="./images/paper/cover2.png" alt="cover2" style="width:450px;"/>
+<br>
+<img src="./images/paper/cover3.png" alt="cover3" style="width:450px;"/>
+<br>
+<img src="./images/paper/cover4.png" alt="cover4" style="width:800px;"/>
+<br>
+<img src="./images/paper/cover5.png" alt="cover5" style="width:800px;"/>
+<br>
+<img src="./images/paper/cover6.png" alt="cover6" style="width:400px;"/>
+<br>
+<img src="./images/paper/cover7.png" alt="cover7" style="width:400px;"/>
+<img src="./images/paper/cover8.png" alt="cover8" style="width:400px;"/>

examples/03.IAXO/BabyIAXO_2024.rml

@@ -0,0 +1,110 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<BabyIAXO>
+	<globals>
+		<constant name="Pi" value="3.1415927" />
+        <constant name="Bores" value="2" />
+        <constant name="SpotRadius" value="0.8" />
+		<constant name="BField" value="1.69" />
+        <constant name="Years" value="1.5" />
+        <constant name="MagnetRadius" value="35" />
+        <constant name="MagnetLength" value="11" />
+		<constant name="GasLength" value="5" />
+        <constant name="WindowEff" value="0.81" />
+        <constant name="OpticsEff" value="0.47*0.9" />
+        <constant name="DetectorEff" value="0.66" />
+        <constant name="BckRate" value="1e-7" />
+	</globals>
+
+	<TRestSensitivity name="VacuumPhase" strategy="nodes" >
+
+		<TRestExperiment name="Vacuum" exposureTime="Years*365*12hr" useAverage="true">
+			<!-- Background -->
+			<TRestComponentFormula name="Flat7" nature="background" >
+				<parameter name="formulaUnits" value="keV^-1" />
+				<cVariable name="energy" range="(0,10)keV" bins="20" />
+				<formula name="bck" expression="Bores*BckRate*Pi*SpotRadius*SpotRadius" />
+			</TRestComponentFormula>
+
+			<!-- Signal -->
+			<TRestAxionHelioscopeSignal name="Vacuum" nature="signal" parameter="ma"
+				conversionType="IAXO" bores="Bores"
+				magnetRadius="MagnetRadius cm" magnetLength="MagnetLength m" magnetStrength="BField T"
+				opticsEfficiency="OpticsEff" windowEfficiency="WindowEff" >
+
+				<cVariable name="energy" range="(0,10)keV" bins="20" />
+
+				<parameter name="firstParameterValue" value="0.001" />
+				<parameter name="lastParameterValue" value="10" />
+				<parameter name="stepParameterValue" value="1.02" />
+				<parameter name="exponential" value="true" />
+
+				<TRestAxionSolarQCDFlux name="LennertHoofPrimakoff" couplingType="g_ag" couplingStrength="1.e-10"
+					fluxDataFile="Primakoff_LennertHoof_202203.dat" seed="137" />
+			</TRestAxionHelioscopeSignal>
+		</TRestExperiment>
+
+	</TRestSensitivity>
+
+	<TRestSensitivity name="CombinedPhase" strategy="nodes" >
+		<TRestExperiment name="Vacuum" exposureTime="1.5*300*12hr" useAverage="true">
+			<!-- Background -->
+			<TRestComponentFormula name="Flat7" nature="background" >
+				<parameter name="formulaUnits" value="keV^-1" />
+				<cVariable name="energy" range="(0,10)keV" bins="20" />
+				<formula name="bck" expression="Bores*BckRate*Pi*SpotRadius*SpotRadius" />
+			</TRestComponentFormula>
+
+			<!-- Signal -->
+			<TRestAxionHelioscopeSignal name="Vacuum" nature="signal" parameter="ma"
+				conversionType="IAXO" bores="Bores"
+				magnetRadius="MagnetRadius cm" magnetLength="MagnetLength m" magnetStrength="BField T"
+				opticsEfficiency="OpticsEff" windowEfficiency="WindowEff" detectorEfficiency="DetectorEff"
+				gasLength="GasLength m" >
+
+				<cVariable name="energy" range="(0,10)keV" bins="20" />
+
+				<TRestAxionSolarQCDFlux name="LennertHoofPrimakoff" couplingType="g_ag" couplingStrength="1.e-10"
+					fluxDataFile="Primakoff_LennertHoof_202203.dat" seed="137" />
+
+
+				<parameter name="firstParameterValue" value="0.001" />
+				<parameter name="lastParameterValue" value="10" />
+				<parameter name="stepParameterValue" value="1.02" />
+				<parameter name="exponential" value="true" />
+			</TRestAxionHelioscopeSignal>
+		</TRestExperiment>
+
+		<TRestExperimentList name="GasPhase" exposureTime="0"
+			componentPattern="output/SignalsBabyIAXO_2024.root" experimentsFile="output/KSVZ.settings" useAverage="true">
+
+			<TRestComponentFormula name="Flat7" nature="background" >
+				<parameter name="formulaUnits" value="keV^-1" />
+				<cVariable name="energy" range="(0,10)keV" bins="20" />
+				<formula name="bck" expression="Bores*BckRate*Pi*SpotRadius*SpotRadius" />
+			</TRestComponentFormula>
+		</TRestExperimentList>
+	</TRestSensitivity>
+
+	<TRestAxionHelioscopeSignal name="GasSignal" nature="signal" parameter="ma"
+		conversionType="IAXO" bores="Bores"
+		magnetRadius="MagnetRadius cm" magnetLength="MagnetLength m" magnetStrength="BField T"
+		opticsEfficiency="OpticsEff" windowEfficiency="WindowEff" detectorEfficiency="DetectorEff"
+		gasLength="GasLength m" >
+
+		<cVariable name="energy" range="(0,10)keV" bins="20" />
+
+		<parameter name="firstParameterValue" value="0.001" />
+		<parameter name="lastParameterValue" value="10" />
+		<parameter name="stepParameterValue" value="1.02" />
+		<parameter name="exponential" value="true" />
+
+		<TRestAxionSolarQCDFlux name="LennertHoofPrimakoff" couplingType="g_ag" couplingStrength="1.e-10"
+			fluxDataFile="Primakoff_LennertHoof_202203.dat" seed="137" />
+
+		<TRestAxionBufferGas name="helium" >
+			<gas name="He" density="0.0025e-6g/cm^3"/>
+		</TRestAxionBufferGas>
+
+	</TRestAxionHelioscopeSignal>
+
+</BabyIAXO>

examples/03.IAXO/GenerateSignalComponents.C

@@ -1,3 +1,5 @@
+#include <filesystem>
+
 Double_t Eo = 0.5;  // keV
 Double_t Ef = 10;   // keV
 

examples/03.IAXO/response.rml

@@ -0,0 +1,5 @@
+<resp>
+        <TRestResponse name="XenonNeon" variable="final_energy">
+            <parameter name="filename" value="XenonNeon_50Pct_1.4bar.N150f" />
+        </TRestResponse>
+</resp>

inc/TRestAxionHelioscopeSignal.h

@@ -53,6 +53,12 @@ class TRestAxionHelioscopeSignal : public TRestComponent {
     /// If an x-ray window is present we may add an efficiency for Ngamma calculation
     Double_t fWindowEfficiency = 1;
 
+    /// The overall detector efficiency
+    Double_t fDetectorEfficiency = 1;
+
+    /// The additional gas length photons need to travel to reach a vacuum region (mm)
+    Double_t fGasLength = 0;
+
     /// It defines the gas mixture we use inside our magnetic field. Vacuum if it is nullptr
     TRestAxionBufferGas* fGas = nullptr;
 
@@ -80,19 +86,21 @@ class TRestAxionHelioscopeSignal : public TRestComponent {
     void SetMagnetRadius(const Double_t& radius) { fMagnetRadius = radius; }
     void SetMagnetLength(const Double_t& length) { fMagnetLength = length; }
     void SetMagnetStrength(const Double_t& strength) { fMagnetStrength = strength; }
+    void SetGasLength(const Double_t& length) { fGasLength = length; }
     void SetType(const std::string& type) { fConversionType = type; }
 
     Int_t GetNumberOfBores() const { return fBores; }
     Double_t GetMagnetRadius() const { return fMagnetRadius; }
     Double_t GetMagnetLength() const { return fMagnetLength; }
     Double_t GetMagnetStrength() const { return fMagnetStrength; }
+    Double_t GetGasLength() const { return fGasLength; }
     std::string GetType() const { return fConversionType; }
 
     void Initialize() override;
     TRestAxionHelioscopeSignal(const char* cfgFileName, const std::string& name);
     TRestAxionHelioscopeSignal();
     ~TRestAxionHelioscopeSignal();
 
-    ClassDefOverride(TRestAxionHelioscopeSignal, 1);
+    ClassDefOverride(TRestAxionHelioscopeSignal, 2);
 };
 #endif

inc/TRestAxionMagneticField.h

@@ -118,6 +118,7 @@ class TRestAxionMagneticField : public TRestMetadata {
     void ReMap(const size_t& n, const TVector3& newMapSize);
 
     void SetTrack(const TVector3& position, const TVector3& direction);
+    void SetUserTrack(const TVector3& start, const TVector3& end);
 
     Double_t GetTrackLength() const { return fTrackLength; }
     TVector3 GetTrackStart() const { return fTrackStart; }

inc/TRestAxionOpticsMirror.h

@@ -99,6 +99,8 @@ class TRestAxionOpticsMirror : public TRestMetadata {
 
     TCanvas* DrawOpticsProperties(std::string options = "", Double_t lowRange = 1.e-9,
                                   Double_t lowRange2 = 1.e-3);
+    TPad* DrawOpticsPropertiesLinear(std::string options = "", Double_t lowRange = 1.e-9,
+                                     Double_t lowRange2 = 1.e-3);
 
     TRestAxionOpticsMirror();
     TRestAxionOpticsMirror(const char* cfgFileName, std::string name = "");

inc/TRestAxionTrueWolterOptics.h

@@ -143,6 +143,13 @@ class TRestAxionTrueWolterOptics : public TRestAxionOptics {
         return 0;
     }
 
+    /// It returns the value of max entrance radius in mm
+    Double_t GetMaxEntranceRadius() { return GetR1().back(); }
+
+    /// It returns the value of min entrance radius in mm
+    Double_t GetMinEntranceRadius() { return GetR1().front(); }
+
+    /// It returns the value min/max entrance radius in mm as a std::pair
     std::pair<Double_t, Double_t> GetRadialLimits() override {
         std::pair<Double_t, Double_t> result(0, 0);
         if (!fR1.empty()) {
@@ -151,6 +158,8 @@ class TRestAxionTrueWolterOptics : public TRestAxionOptics {
         return result;
     }
 
+    TRestSpiderMask* const& GetSpiderMask() const { return fSpiderMask; }
+
     void SetMirror() override {
         Double_t x = fEntrancePosition.X();
         Double_t y = fEntrancePosition.Y();

macros/REST_Axion_AccurateEfficiencies.C

@@ -0,0 +1,175 @@
+#include "TCanvas.h"
+#include "TGraph.h"
+#include "TLatex.h"
+#include "TLegend.h"
+#include "TLine.h"
+#include "TRestAxionBufferGas.h"
+#include "TRestAxionField.h"
+
+Double_t fromEnergy = 0.5;
+Double_t toEnergy = 10;
+
+Double_t incidenceAngle = 0.2;
+Double_t deltaE = 0.1;
+
+//*******************************************************************************************************
+//*** Description: This script will launch the integration of the axion-field with given parameters.
+//*** It allows to test different magnetic field cell sizes, for a given mass that can be off-resonance
+//*** for dm different from zero, and a given maximum tolerance or error for the integration routine.
+//***
+//*** The macro sets the TRestAxionField under debug mode to print the different results on screen.
+//***
+//*** --------------
+//*** Usage: restManager FieldIntegrationTests [sX=10] [sX=10] [sZ=10] [dm=0.01] [tolerance=0.1] [Ea=4.2]
+//*** --------------
+//***
+//*** Author: Javier Galan
+//*******************************************************************************************************
+int REST_Axion_AccurateEfficiencies(std::string fluxFile = "fluxes.rml",
+                                    std::string fluxName = "LennertHoofPrimakoff",
+                                    std::string opticsFile = "xmmTrueWolter.rml",
+                                    std::string opticsName = "xmm") {
+    TRestAxionTrueWolterOptics optics(opticsFile.c_str(), opticsName.c_str());
+    TRestAxionOpticsMirror* mirror = optics.GetMirrorProperties();
+
+    TRestAxionSolarQCDFlux flux(fluxFile.c_str(), fluxName.c_str());
+    flux.Initialize();
+
+    Double_t R2sum = 0;
+    Double_t fluxSum = 0;
+    Double_t maxFlux = 0;
+    for (Double_t energy = fromEnergy; energy < toEnergy; energy += deltaE) {
+        Double_t R = mirror->GetReflectivity(incidenceAngle, energy);
+        fluxSum += flux.GetFluxAtEnergy(energy, 0);
+        if (flux.GetFluxAtEnergy(energy, 0) > maxFlux) maxFlux = flux.GetFluxAtEnergy(energy, 0);
+        R2sum += flux.GetFluxAtEnergy(energy, 0) * R * R;
+    }
+
+    Double_t R2eff = R2sum / fluxSum;
+    std::cout << "R2eff: " << R2eff << std::endl;
+
+    Double_t Rout = optics.GetMaxEntranceRadius();
+    Double_t Rin = optics.GetMinEntranceRadius();
+
+    TRestSpiderMask* sMask = optics.GetSpiderMask();
+    Double_t na = sMask->GetNumberOfArms();
+    Double_t wa = sMask->GetArmsWidth();
+
+    std::vector<Double_t> r = optics.GetR1();
+    std::vector<Double_t> th = optics.GetThickness();
+
+    Double_t Aeff = (TMath::Pi() - 0.5 * na * wa) * (Rout * Rout - Rin * Rin);
+    for (size_t n = 0; n < r.size(); n++) Aeff -= (2 * TMath::Pi() - na * wa) * r[n] * th[n];
+
+    std::cout << "Aeff (optics): " << Aeff / Rout / Rout / TMath::Pi() << std::endl;
+
+    TCanvas c;
+    c.SetCanvasSize(2400, 1800);
+    c.SetWindowSize(2400, 1800);
+    c.Divide(2, 1);
+
+    c.cd(1);
+    optics.GetMirrorProperties()->DrawOpticsPropertiesLinear();
+
+    c.cd(2);
+    optics.DrawParticleTracks();
+
+    c.Print("optics.pdf");
+
+    /// Extracted from x-ray window MicromegasStrongBack
+    na = 8;
+    wa = 2.64 * TMath::Pi() / 180.;
+    Rout = 8.5;
+    Rin = 4.55;
+    Double_t Ro = 4.25;
+
+    Aeff = (TMath::Pi() - 0.5 * na * wa) * (Rout * Rout - Rin * Rin) + TMath::Pi() * Ro * Ro;
+
+    std::cout << "Aeff (window): " << Aeff / Rout / Rout / TMath::Pi() << std::endl;
+
+    TRestAxionXrayWindow strongBack("windows.rml", "MicromegasStrongBack");
+    TRestAxionXrayWindow mylar("windows.rml", "MicromegasMylar");
+    TRestAxionXrayWindow aluminum("windows.rml", "MicromegasAluminumFoil");
+
+    TGraph* mylarGraph = new TGraph();
+    mylarGraph->SetName("Mylar");
+    mylarGraph->SetLineColor(49);
+    mylarGraph->SetLineWidth(2);
+
+    TGraph* aluminumGraph = new TGraph();
+    aluminumGraph->SetName("Aluminum");
+    aluminumGraph->SetLineColor(46);
+    aluminumGraph->SetLineWidth(2);
+
+    TGraph* solarGraph = new TGraph();
+    solarGraph->SetName("SolarFlux");
+    solarGraph->SetLineColor(43);
+    solarGraph->SetLineWidth(2);
+
+    Double_t WeffSum = 0;
+    Double_t AlSum = 0;
+    Double_t MySum = 0;
+    for (Double_t energy = deltaE; energy < toEnergy; energy += deltaE) {
+        Double_t tMy = mylar.GetTransmission(energy, 0, 0);
+        Double_t tAl = aluminum.GetTransmission(energy, 0, 0);
+
+        mylarGraph->SetPoint(mylarGraph->GetN(), energy, tMy);
+        aluminumGraph->SetPoint(aluminumGraph->GetN(), energy, tAl);
+        solarGraph->SetPoint(solarGraph->GetN(), energy, flux.GetFluxAtEnergy(energy, 0) / maxFlux);
+
+        WeffSum += flux.GetFluxAtEnergy(energy, 0) * tMy * tAl;
+        AlSum += flux.GetFluxAtEnergy(energy, 0) * tAl;
+        MySum += flux.GetFluxAtEnergy(energy, 0) * tMy;
+    }
+
+    WeffSum = WeffSum / fluxSum;
+    AlSum = AlSum / fluxSum;
+    MySum = MySum / fluxSum;
+    std::cout << "AlSum: " << AlSum << std::endl;
+    std::cout << "MySum: " << MySum << std::endl;
+    std::cout << "WeffSum: " << WeffSum << std::endl;
+
+    TCanvas c2;
+    c2.SetCanvasSize(1200, 900);
+    c2.SetWindowSize(1200, 900);
+    // c2.SetLogy();
+
+    TPad* pad2 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97);
+    //   pad1->Divide(2, 2);
+    pad2->SetLogy();
+    pad2->Draw();
+
+    ////// Drawing reflectivity versus angle
+    pad2->SetRightMargin(0.09);
+    pad2->SetLeftMargin(0.25);
+    pad2->SetBottomMargin(0.15);
+
+    mylarGraph->GetXaxis()->SetLimits(0, 10);
+    //   mylarGraph->GetHistogram()->SetMaximum(1);
+    mylarGraph->GetHistogram()->SetMinimum(0);
+
+    mylarGraph->GetXaxis()->SetTitle("Energy [keV]");
+    mylarGraph->GetXaxis()->SetTitleSize(0.04);
+    mylarGraph->GetXaxis()->SetLabelSize(0.04);
+    mylarGraph->GetYaxis()->SetTitle("Transmission");
+    mylarGraph->GetYaxis()->SetTitleOffset(1.2);
+    mylarGraph->GetYaxis()->SetTitleSize(0.04);
+    mylarGraph->GetYaxis()->SetLabelSize(0.04);
+    mylarGraph->Draw("AL");
+    aluminumGraph->Draw("L");
+    // solarGraph->Draw("L");
+
+    Double_t lx1 = 0.6, ly1 = 0.55, lx2 = 0.8, ly2 = 0.75;
+    TLegend* legend = new TLegend(lx1, ly1, lx2, ly2);
+    legend->SetTextSize(0.03);
+    //   legend->SetHeader("Widnows", "C");  // option "C" allows to center the header
+    legend->AddEntry("Mylar", "Mylar", "l");
+    legend->AddEntry("Aluminum", "Aluminum", "l");
+    legend->Draw();
+
+    c2.Print("windows.pdf");
+
+    c.Print("tracks.png");
+
+    return 0;
+}