TRestAxionTransmissionProcess. Fixing the location of the window

.github/workflows/validation.yml

@@ -43,11 +43,7 @@ jobs:
         with:
           cmake-flags: "-DCMAKE_INSTALL_PREFIX=${{ env.REST_PATH }} -DCMAKE_BUILD_TYPE=${{ env.CMAKE_BUILD_TYPE }} -DREST_WELCOME=ON -DRESTLIB_AXION=ON"
           branch: ${{ env.BRANCH_NAME }}
-      - name: Load REST libraries
-        run: |
-          source ${{ env.REST_PATH }}/thisREST.sh
-          restRoot -b -q
-      - name: Cache axionlib installation
+      - name: Cache framework installation
         id: axionlib-install-cache
         uses: actions/cache@v3
         with:
@@ -144,3 +140,14 @@ jobs:
         wget https://sultan.unizar.es/axionlib-data/opticsMirror/Transmission_Single_Au_250_Ni_0.4.N901f
         restManager --c opticsBench.rml
         restRoot ValidateXMM.C'("OpticsBench_Yaw_0.02_Dev_0.005_BabyIAXO_Run00100.root")'
+    - name: Window transmission
+      run: |
+        source ${{ env.REST_PATH }}/thisREST.sh
+        cd pipeline/ray-tracing/transmission/
+        wget https://sultan.unizar.es/axionlib-data/transmission/Si.sol
+        restManager --c emptyWindow.rml
+        restRoot ValidateTransmission.C'("EmptyVacuumWindow.root", 0.07, 0.09)'
+        restManager --c shiftedWindow.rml
+        restRoot ValidateTransmission.C'("ShiftedVacuumWindow", 0.3, 0.5)'
+        restManager --c siWindow.rml
+        restRoot ValidateTransmission.C'("SiWindow.root", 0.9, 1)'

README.md

@@ -3,6 +3,7 @@
 [![website](https://img.shields.io/badge/user-guide-E8B6FF.svg)](https://rest-for-physics.github.io)
 [![api](https://img.shields.io/badge/user-API-FFCA78.svg)](https://sultan.unizar.es/rest/)
 [![forum](https://img.shields.io/badge/user-forum-AAFF90.svg)](https://rest-forum.unizar.es/)
+[![github pipeline](https://github.com/rest-for-physics/axionlib/actions/workflows/validation.yml/badge.svg)](https://github.com/rest-for-physics/axionlib/commits/master)
 
 This is a REST-for-Physics library used to generate solar axions and obtain the detection probability function. The library allows to define a generic helioscope setup - buffer gas, magnetic field, optics response, photon transmission - through REST metadata structures.
 

inc/TRestAxionXrayWindow.h

@@ -29,14 +29,11 @@
 //! A metadata class to create x-ray transmission window definitions
 class TRestAxionXrayWindow : public TRestMetadata {
    private:
-    /// Position of the center of the window in mm
-    TVector3 fCenter = TVector3(0, 0, 0);  //<
-
     /// Thicknesss of window material in mm
     Double_t fThickness = 0.0;  //<
 
     /// Window material name
-    std::string fMaterial = "Si";  //<
+    std::string fMaterial = "Vacuum";  //<
 
     /// A mask defining a pattern where the transmission will be effective
     TRestPatternMask* fMask = nullptr;  //<
@@ -75,6 +72,6 @@ class TRestAxionXrayWindow : public TRestMetadata {
 
     ~TRestAxionXrayWindow();
 
-    ClassDef(TRestAxionXrayWindow, 1);
+    ClassDef(TRestAxionXrayWindow, 2);
 };
 #endif

pipeline/metadata/transmission/windowPlot.py

@@ -27,8 +27,11 @@
 
 totalSamples = 100000
 
+print ("Loading cathode" )
 cathode = ROOT.TRestAxionXrayWindow("windows.rml", "cathode")
+print ("Loading strongBack" )
 strongBack = ROOT.TRestAxionXrayWindow("windows.rml", patternType)
+print ("Loading silicon foil" )
 siFoil = ROOT.TRestAxionXrayWindow("windows.rml", "siliconFoil")
 
 radius = strongBack.GetWindowRadius()
@@ -136,6 +139,9 @@
 
 c1.Print(outfname)
 
+print ( "Low: " + str( histL.Integral() ) )
+print ( "Mid: " + str( histM.Integral() ) )
+print ( "High: " + str( histH.Integral() ) )
 if( histL.Integral() < 10000 ):
     print( "Effective counts at low energy below 10000!!")
     print ( "Low: " + str( histL.Integral() ) )

pipeline/ray-tracing/transmission/ValidateTransmission.C

@@ -0,0 +1,34 @@
+#include <TH1D.h>
+#include <TRestRun.h>
+
+Int_t ValidateTransmission(std::string fname, Double_t lowEff = 0.07, Double_t highEff = 0.09) {
+    std::cout << "Filename : " << fname << std::endl;
+    TRestRun* run = new TRestRun(fname);
+
+    if (run->GetEntries() != 10000) {
+        std::cout << "Error. Number of entries is not 10000!" << std::endl;
+
+        return 1;
+    }
+
+    run->GetAnalysisTree()->Draw("bore_transmission", "bore_transmission");
+
+    TH1D* h = (TH1D*)run->GetAnalysisTree()->GetHistogram();
+
+    if (h == nullptr) {
+        std::cout << "Problems generating histogram" << std::endl;
+        return 2;
+    }
+
+    Double_t integral = h->Integral();
+    std::cout << "Photon detection efficiency : " << integral / run->GetEntries() << std::endl;
+    Double_t efficiency = integral / run->GetEntries();
+    if (efficiency < lowEff || efficiency > highEff) {
+        std::cout << "Wrong number of photons detected!" << std::endl;
+        return 3;
+    }
+
+    delete run;
+
+    return 0;
+}

pipeline/ray-tracing/transmission/emptyWindow.rml

@@ -0,0 +1,32 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<TRestManager name="RESTManagerSim" title="Template manager to process a simulation generated by restG4." verboseLevel="info">
+    <globals>
+        <parameter name="verboseLevel" value="essential"/>
+    </globals>
+    <TRestRun name="axionRun" title="BabyIAXO V1.0" verboseLevel="info">
+        <parameter name="experimentName" value="BabyIAXO"/>
+        <parameter name="runType" value="Transmission"/>
+        <parameter name="runNumber" value="103"/>
+        <parameter name="runTag" value="emptyWindow"/>
+        <parameter name="runDescription" value=""/>
+        <parameter name="user" value="${USER}"/>
+        <parameter name="verboseLevel" value="2"/>
+        <parameter name="outputFileName" value="EmptyVacuumWindow.root"/>
+        <TRestAxionXrayWindow name="magnetBore" material="Vacuum">
+            <TRestPatternMask name="mask" maskRadius="10cm" verboseLevel="info"/>
+        </TRestAxionXrayWindow>
+    </TRestRun>
+    <TRestProcessRunner name="EventProcess" verboseLevel="info">
+        <parameter name="eventsToProcess" value="10000"/>
+        <addProcess type="TRestAxionGeneratorProcess" name="axionGen">
+            <parameter name="generatorType" value="flat"/>
+            <parameter name="maxEnergy" value="10"/>
+            <parameter name="targetRadius" value="35cm"/>
+        </addProcess>
+        <addProcess type="TRestAxionTransmissionProcess" name="bore" center="(0,0,0)">
+            <window name="magnetBore"/>
+        </addProcess>
+        <addProcess type="TRestAxionAnalysisProcess" name="axGen1" observables="all"/>
+    </TRestProcessRunner>
+    <addTask command="EventProcess-&gt;RunProcess()" value="ON"/>
+</TRestManager>

pipeline/ray-tracing/transmission/shiftedWindow.rml

@@ -0,0 +1,32 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<TRestManager name="RESTManagerSim" title="Template manager to process a simulation generated by restG4." verboseLevel="info">
+    <globals>
+        <parameter name="verboseLevel" value="essential"/>
+    </globals>
+    <TRestRun name="axionRun" title="BabyIAXO V1.0" verboseLevel="info">
+        <parameter name="experimentName" value="BabyIAXO"/>
+        <parameter name="runType" value="Transmission"/>
+        <parameter name="runNumber" value="103"/>
+        <parameter name="runTag" value="shiftedWindow"/>
+        <parameter name="runDescription" value=""/>
+        <parameter name="user" value="${USER}"/>
+        <parameter name="verboseLevel" value="2"/>
+        <parameter name="outputFileName" value="ShiftedVacuumWindow.root"/>
+        <TRestAxionXrayWindow name="magnetBore" material="Vacuum">
+            <TRestPatternMask name="mask" maskRadius="35cm" verboseLevel="info"/>
+        </TRestAxionXrayWindow>
+    </TRestRun>
+    <TRestProcessRunner name="EventProcess" verboseLevel="info">
+        <parameter name="eventsToProcess" value="10000"/>
+        <addProcess type="TRestAxionGeneratorProcess" name="axionGen">
+            <parameter name="generatorType" value="flat"/>
+            <parameter name="maxEnergy" value="10"/>
+            <parameter name="targetRadius" value="35cm"/>
+        </addProcess>
+        <addProcess type="TRestAxionTransmissionProcess" name="bore" center="(35,0,0)cm">
+            <window name="magnetBore"/>
+        </addProcess>
+        <addProcess type="TRestAxionAnalysisProcess" name="axGen1" observables="all"/>
+    </TRestProcessRunner>
+    <addTask command="EventProcess-&gt;RunProcess()" value="ON"/>
+</TRestManager>

pipeline/ray-tracing/transmission/siWindow.rml

@@ -0,0 +1,32 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<TRestManager name="RESTManagerSim" title="Template manager to process a simulation generated by restG4." verboseLevel="info">
+    <globals>
+        <parameter name="verboseLevel" value="essential"/>
+    </globals>
+    <TRestRun name="axionRun" title="BabyIAXO V1.0" verboseLevel="info">
+        <parameter name="experimentName" value="BabyIAXO"/>
+        <parameter name="runType" value="Transmission"/>
+        <parameter name="runNumber" value="103"/>
+        <parameter name="runTag" value="SiVacuumWindow"/>
+        <parameter name="runDescription" value=""/>
+        <parameter name="user" value="${USER}"/>
+        <parameter name="verboseLevel" value="2"/>
+        <parameter name="outputFileName" value="SiWindow.root"/>
+        <TRestAxionXrayWindow name="siWindow" material="Si" thickness="50nm">
+            <TRestPatternMask name="mask" maskRadius="35cm" verboseLevel="info"/>
+        </TRestAxionXrayWindow>
+    </TRestRun>
+    <TRestProcessRunner name="EventProcess" verboseLevel="info">
+        <parameter name="eventsToProcess" value="10000"/>
+        <addProcess type="TRestAxionGeneratorProcess" name="axionGen">
+            <parameter name="generatorType" value="flat"/>
+            <parameter name="maxEnergy" value="10"/>
+            <parameter name="targetRadius" value="35cm"/>
+        </addProcess>
+        <addProcess type="TRestAxionTransmissionProcess" name="bore" center="(0,0,0)">
+            <window name="siWindow"/>
+        </addProcess>
+        <addProcess type="TRestAxionAnalysisProcess" name="axGen1" observables="all"/>
+    </TRestProcessRunner>
+    <addTask command="EventProcess-&gt;RunProcess()" value="ON"/>
+</TRestManager>

src/TRestAxionTransmissionProcess.cxx

@@ -110,6 +110,8 @@ void TRestAxionTransmissionProcess::Initialize() {
 void TRestAxionTransmissionProcess::InitProcess() {
     RESTDebug << "Entering ... TRestAxionGeneratorProcess::InitProcess" << RESTendl;
 
+    RegisterAllObservables();
+
     fXrayWindows.clear();
     for (const auto& wName : fWindowNames) {
         TRestAxionXrayWindow* w = (TRestAxionXrayWindow*)GetMetadata(wName);
@@ -133,11 +135,16 @@ TRestEvent* TRestAxionTransmissionProcess::ProcessEvent(TRestEvent* evInput) {
     fAxionEvent = (TRestAxionEvent*)evInput;
 
     TVector3 inPos = fAxionEvent->GetPosition();
+    TVector3 inDir = fAxionEvent->GetDirection();
+
+    /// The component is placed at (0,0,0). It is TRestAxionEventProcess the responsible to translate the
+    /// component (in reality the particle) according to fCenter.
+    TVector3 newPos = REST_Physics::MoveToPlane(inPos, inDir, TVector3(0, 0, 1), TVector3(0, 0, 0));
 
     Double_t transmission = 1;
-    Double_t x = inPos.X();
-    Double_t y = inPos.Y();
-    Double_t z = inPos.Z();
+    Double_t x = newPos.X();
+    Double_t y = newPos.Y();
+    Double_t z = newPos.Z();
     Double_t en = fAxionEvent->GetEnergy();
 
     RESTDebug << "Particle position to evaluate window transmission. " << RESTendl;
@@ -146,6 +153,7 @@ TRestEvent* TRestAxionTransmissionProcess::ProcessEvent(TRestEvent* evInput) {
     for (const auto& window : fXrayWindows) {
         transmission *= window->GetTransmission(en, x, y);
     }
+    RESTDebug << "Transmission: " << transmission << RESTendl;
 
     SetObservableValue("transmission", transmission);
 

src/TRestAxionXrayWindow.cxx

@@ -196,14 +196,14 @@ void TRestAxionXrayWindow::ReadMaterial() {
 /// not hitted.
 ///
 Double_t TRestAxionXrayWindow::GetTransmission(Double_t energy, Double_t x, Double_t y) {
-    if (fEnergy.size() == 0) ReadMaterial();
+    if (fMaterial != "Vacuum" && fEnergy.size() == 0) ReadMaterial();
 
-    Double_t xNew = x - fCenter.X();
-    Double_t yNew = y - fCenter.Y();
+    if (fMask && x * x + y * y > fMask->GetMaskRadius() * fMask->GetMaskRadius()) return 0;
 
-    if (fMask && xNew * xNew + yNew * yNew > fMask->GetMaskRadius() * fMask->GetMaskRadius()) return 0;
+    if (fMask && !fMask->HitsPattern(x, y)) return 1.;
 
-    if (fMask && !fMask->HitsPattern(xNew, yNew)) 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);
 
@@ -310,7 +310,12 @@ void TRestAxionXrayWindow::InitFromConfigFile() {
         delete fMask;
         fMask = nullptr;
     }
-    fMask = (TRestPatternMask*)this->InstantiateChildMetadata("Mask", "mask");
+    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;
+    }
 }
 
 ///////////////////////////////////////////////
@@ -319,8 +324,6 @@ void TRestAxionXrayWindow::InitFromConfigFile() {
 void TRestAxionXrayWindow::PrintMetadata() {
     TRestMetadata::PrintMetadata();
 
-    RESTMetadata << "Window center: ( " << fCenter.X() << ", " << fCenter.Y() << ", " << fCenter.Z() << ") mm"
-                 << RESTendl;
     RESTMetadata << "Thickness: " << fThickness * units("um") << " um" << RESTendl;
     RESTMetadata << "Material: " << fMaterial << RESTendl;
     RESTMetadata << "----" << RESTendl;