Merge pull request #93 from rest-for-physics/jgalan_bIAXO_example

Helioscope signal sensitivity
rest-for-physics · Jun 11, 2024 · f6cfc9a · f6cfc9a
2 parents 3802d10 + 79e2cbd
commit f6cfc9a
Show file tree

Hide file tree

Showing 35 changed files with 4,004 additions and 16 deletions.
diff --git a/.github/workflows/validation.yml b/.github/workflows/validation.yml
@@ -266,3 +266,27 @@ jobs:
         export REST_GAS_STEP=2
         restManager --c BabyIAXO_Micromegas_XMM.rml --o output.root
         restRoot -q ValidateChain.C'("output.root")'
+
+  Sensitivity:
+    name: Sensitivity tests
+    runs-on: ubuntu-latest
+    container:
+      image: ghcr.io/lobis/root-geant4-garfield:rest-for-physics
+    needs: [ build-axionlib ]
+    steps:
+    - uses: rest-for-physics/framework/.github/actions/checkout@master
+      with:
+        branch: ${{ env.BRANCH_NAME }}
+        repository: rest-for-physics/axionlib
+        path: ${{ env.AXION_LIB_PATH }}
+    - name: Restore cache
+      uses: actions/cache@v3
+      id: axionlib-install-cache
+      with:
+        key: ${{ env.BRANCH_NAME }}-${{ github.sha }}
+        path: ${{ env.REST_PATH }}
+    - name: Dummy IAXO
+      run: |
+        source ${{ env.REST_PATH }}/thisREST.sh
+        cd ${{ env.AXION_LIB_PATH }}/pipeline/sensitivity/
+        restRoot -b -q Validation.C
diff --git a/examples/full-ray-tracing/helioscope.rml → examples/01.full-ray-tracing/helioscope.rml b/examples/full-ray-tracing/helioscope.rml → examples/01.full-ray-tracing/helioscope.rml
diff --git a/examples/full-ray-tracing/summary.rml → examples/01.full-ray-tracing/summary.rml b/examples/full-ray-tracing/summary.rml → examples/01.full-ray-tracing/summary.rml
diff --git a/examples/optics/opticsBench.rml → examples/02.optics/opticsBench.rml b/examples/optics/opticsBench.rml → examples/02.optics/opticsBench.rml
diff --git a/examples/optics/opticsBenchTrueW1.rml → examples/02.optics/opticsBenchTrueW1.rml b/examples/optics/opticsBenchTrueW1.rml → examples/02.optics/opticsBenchTrueW1.rml
diff --git a/examples/optics/opticsPlots.rml → examples/02.optics/opticsPlots.rml b/examples/optics/opticsPlots.rml → examples/02.optics/opticsPlots.rml
diff --git a/examples/optics/opticsPlotsAllTransm.rml → examples/02.optics/opticsPlotsAllTransm.rml b/examples/optics/opticsPlotsAllTransm.rml → examples/02.optics/opticsPlotsAllTransm.rml
diff --git a/examples/03.IAXO/BabyIAXO.rml b/examples/03.IAXO/BabyIAXO.rml
@@ -0,0 +1,106 @@
+<?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.309" />
+		<constant name="BField" value="1.73" />
+		<constant name="Years" value="1.5" />
+		<constant name="MagnetRadius" value="35" />
+		<constant name="MagnetLength" value="10" />
+		<constant name="WindowEff" value="0.7" />
+		<constant name="OpticsEff" value="0.35" />
+		<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" >
+
+				<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.root" experimentsFile="output/BabyIAXO.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" >
+
+		<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>
diff --git a/examples/03.IAXO/GenerateSignalComponents.C b/examples/03.IAXO/GenerateSignalComponents.C
@@ -0,0 +1,64 @@
+Double_t Eo = 0.5;  // keV
+Double_t Ef = 10;   // keV
+
+Int_t GenerateSignalComponents(std::string rmlFile, std::string name, Double_t totalYears = 1.5,
+                               size_t firstSteps = 5, Double_t firstYears = 0.5) {
+    Double_t totalExposureTime = totalYears * 365. * 12 * 3600;
+    Double_t firstStepsExposure = firstYears * 365. * 12 * 3600;
+
+    TRestAxionField field;
+    std::vector<std::pair<Double_t, Double_t>> scanSteps =
+        field.GetMassDensityScanning("He", 0.25, 20);  // Up to 0.25 eV
+
+    TRestAxionHelioscopeSignal gasPhase(rmlFile.c_str(), name.c_str());
+    std::filesystem::path filePath = rmlFile;
+    std::filesystem::path newExtension = ".settings";
+    std::string settingsFile = filePath.replace_extension(newExtension);
+    std::filesystem::path rootExtension = ".root";
+    std::string componentsFile = "Signals" + (std::string)filePath.replace_extension(rootExtension);
+
+    std::filesystem::create_directories("output");
+    settingsFile = "output/" + settingsFile;
+    componentsFile = "output/" + componentsFile;
+
+    TFile* f = TFile::Open((TString)componentsFile, "RECREATE");
+    std::vector<Double_t> exposureTimes;
+    for (size_t n = 0; n < scanSteps.size(); n++) {
+        Double_t mass = scanSteps[n].first;
+        Double_t density = scanSteps[n].second;
+
+        gasPhase.SetName("P" + (TString)IntegerToString(n + 1));
+        gasPhase.GetGas()->SetGasDensity("He", density);
+        gasPhase.RegenerateHistograms();
+
+        Double_t nGamma = gasPhase.GetSignalRate(mass, Eo, Ef);
+        Double_t ksvzFactor = 3.75523 * mass;
+        Double_t exposureTime = 0;
+        if (n < firstSteps)
+            exposureTime = firstStepsExposure / firstSteps;
+        else
+            exposureTime = TMath::Log(20.) / TMath::Power(ksvzFactor, 4) / nGamma;
+        exposureTimes.push_back(exposureTime);
+
+        gasPhase.Write("P" + (TString)IntegerToString(n + 1));
+    }
+    f->Close();
+
+    std::cout << "Generated signals file: " << componentsFile << std::endl;
+
+    Double_t generatedExposureTime = 0;
+    for (size_t n = firstSteps; n < scanSteps.size(); n++) generatedExposureTime += exposureTimes[n];
+
+    FILE* g = fopen(settingsFile.c_str(), "wt");
+    for (size_t n = 0; n < firstSteps; n++) fprintf(g, "%lf\tP%d\n", exposureTimes[n], (int)(n + 1));
+
+    for (size_t n = firstSteps; n < scanSteps.size(); n++)
+        fprintf(g, "%lf\tP%d\n",
+                (totalExposureTime - firstStepsExposure) * exposureTimes[n] / generatedExposureTime,
+                (int)(n + 1));
+    fclose(g);
+
+    std::cout << "Generated settings file: " << settingsFile << std::endl;
+
+    return 0;
+}
diff --git a/examples/03.IAXO/IAXO.rml b/examples/03.IAXO/IAXO.rml
@@ -0,0 +1,105 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<BabyIAXO>
+	<globals>
+		<constant name="Pi" value="3.1415927" />
+		<constant name="Bores" value="8" />
+		<constant name="SpotRadius" value="0.2185" /> <!-- Area 0.15 cm2 -->
+		<constant name="BField" value="2.55" />
+		<constant name="Years" value="3" />
+		<constant name="MagnetRadius" value="30" />
+		<constant name="MagnetLength" value="20" />
+		<constant name="WindowEff" value="0.8" />
+		<constant name="OpticsEff" value="0.7" />
+		<constant name="BckRate" value="1e-8" />
+	</globals>
+
+	<TRestSensitivity name="VacuumPhase" strategy="nodes" >
+
+		<TRestExperiment name="Vacuum" exposureTime="Years*365*12hr" useAverage="true">
+			<!-- Background -->
+			<TRestComponentFormula name="Flat8" 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="Flat8" 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" />
+
+				<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/SignalsIAXO.root" experimentsFile="output/IAXO.settings" useAverage="true">
+
+			<TRestComponentFormula name="Flat8" 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" >
+
+		<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>