diff --git a/.github/workflows/validation.yml b/.github/workflows/validation.yml index 6cdef434..91fb87da 100644 --- a/.github/workflows/validation.yml +++ b/.github/workflows/validation.yml @@ -41,7 +41,7 @@ jobs: - name: Build and install uses: rest-for-physics/framework/.github/actions/build@master with: - cmake-flags: "-DCMAKE_INSTALL_PREFIX=${{ env.REST_PATH }} -DCMAKE_BUILD_TYPE=${{ env.CMAKE_BUILD_TYPE }} -DREST_WELCOME=ON -DRESTLIB_AXION=ON" + cmake-flags: "-DCMAKE_INSTALL_PREFIX=${{ env.REST_PATH }} -DCMAKE_BUILD_TYPE=${{ env.CMAKE_BUILD_TYPE }} -DREST_WELCOME=ON -DRESTLIB_AXION=ON -DREST_MPFR=ON" branch: ${{ env.BRANCH_NAME }} - name: Cache framework installation id: axionlib-install-cache @@ -50,6 +50,55 @@ jobs: key: ${{ env.BRANCH_NAME }}-${{ github.sha }} path: ${{ env.REST_PATH }} + macros-health: + name: "Macros with clean error output" + runs-on: ubuntu-latest + container: + image: ghcr.io/lobis/root-geant4-garfield:rest-for-physics + needs: [ build-axionlib ] + steps: + - uses: actions/checkout@v3 + - name: Restore cache + uses: actions/cache@v3 + id: axionlib-install-cache + with: + key: ${{ env.BRANCH_NAME }}-${{ github.sha }} + path: ${{ env.REST_PATH }} + - name: Load restRootMacros in axiolib standalone + run: | + source ${{ env.REST_PATH }}/thisREST.sh + cd pipeline + export DISPLAY=localhost:0.0 + echo "Running validation script" + python3 validateMacros.py + + Physics: + name: Check physics + runs-on: ubuntu-latest + container: + image: ghcr.io/lobis/root-geant4-garfield:rest-for-physics + needs: [ build-axionlib ] + steps: + - uses: actions/checkout@v3 + - name: Restore cache + uses: actions/cache@v3 + id: axionlib-install-cache + with: + key: ${{ env.BRANCH_NAME }}-${{ github.sha }} + path: ${{ env.REST_PATH }} + - name: Basic physics tests + run: | + source ${{ env.REST_PATH }}/thisREST.sh + cd pipeline/physics/ + wget https://sultan.unizar.es/axionlib-data/bufferGas/He.abs + wget https://sultan.unizar.es/axionlib-data/bufferGas/He.nff + wget https://sultan.unizar.es/axionlib-data/bufferGas/Ne.abs + wget https://sultan.unizar.es/axionlib-data/bufferGas/Ne.nff + wget https://sultan.unizar.es/axionlib-data/bufferGas/bufferGases.rml + wget https://sultan.unizar.es/axionlib-data/magneticField/fields.rml + wget https://sultan.unizar.es/axionlib-data/magneticField/Bykovskiy_201906.dat + restRoot -b -q AxionPhysicsValidation.C + Metadata: name: Check metadata runs-on: ubuntu-latest @@ -134,12 +183,13 @@ jobs: source ${{ env.REST_PATH }}/thisREST.sh export REST_NEVENTS=1000 export REST_RUN=100 + cd pipeline/ray-tracing/optics/ wget https://rest-for-physics.github.io/axionlib-data/optics/xmm.rml wget https://rest-for-physics.github.io/axionlib-data/optics/XMM.Wolter wget https://rest-for-physics.github.io/axionlib-data/opticsMirror/Reflectivity_Single_Au_250_Ni_0.4.N901f wget https://rest-for-physics.github.io/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")' + restRoot -b -q ValidateXMM.C'("OpticsBench_Yaw_0.05_Dev_0.005_BabyIAXO_Run00100.root")' - name: Window transmission run: | source ${{ env.REST_PATH }}/thisREST.sh @@ -150,4 +200,15 @@ jobs: restManager --c shiftedWindow.rml restRoot ValidateTransmission.C'("ShiftedVacuumWindow", 0.3, 0.5)' restManager --c siWindow.rml - restRoot ValidateTransmission.C'("SiWindow.root", 0.9, 1)' + restRoot -b -q ValidateTransmission.C'("SiWindow.root", 0.9, 1)' + - name: Axion-field integration + run: | + source ${{ env.REST_PATH }}/thisREST.sh + cd pipeline/ray-tracing/axion-field/ + wget https://rest-for-physics.github.io/axionlib-data/magneticField/fields.rml + wget https://rest-for-physics.github.io/axionlib-data/magneticField/Bykovskiy_201906.dat + wget https://sultan.unizar.es/axionlib-data/bufferGas/He.abs + wget https://sultan.unizar.es/axionlib-data/bufferGas/He.nff + wget https://rest-for-physics.github.io/axionlib-data/bufferGas/bufferGases.rml + restManager --c photonConversion.rml + restRoot -b -q Validate.C diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index fce91ddb..2e0f6b9c 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -45,7 +45,7 @@ build: - cd ../../../ - mkdir build - cd build - - cmake ../ -DREST_SOLAXFLUX=OFF -DREST_WELCOME=ON -DRESTLIB_AXION=ON -DREST_GARFIELD=OFF -DREST_G4=OFF -DCMAKE_INSTALL_PREFIX=${CI_PROJECT_DIR}/install -DMPFR_PATH=${CI_PROJECT_DIR}/mpfr-4.0.2/install + - cmake ../ -DREST_SOLAXFLUX=OFF -DREST_WELCOME=ON -DRESTLIB_AXION=ON -DREST_GARFIELD=OFF -DREST_G4=OFF -DCMAKE_INSTALL_PREFIX=${CI_PROJECT_DIR}/install -DMPFR_PATH=${CI_PROJECT_DIR}/mpfr-4.0.2/install -DREST_MPFR=ON - make install -j2 # - . ${CI_PROJECT_DIR}/framework/source/libraries/axion/external/solarAxionFlux/bin/thisSolarAxionFluxLib.sh - . ${CI_PROJECT_DIR}/install/thisREST.sh diff --git a/data b/data index d1e52b19..c460ad98 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit d1e52b1968a08b3cb5ea370283e4648bc5da107a +Subproject commit c460ad9853e02c7a15e21f31bf77610f959ab50c diff --git a/examples/full-ray-tracing/helioscope.rml b/examples/full-ray-tracing/helioscope.rml new file mode 100644 index 00000000..5246b30d --- /dev/null +++ b/examples/full-ray-tracing/helioscope.rml @@ -0,0 +1,85 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/full-ray-tracing/summary.rml b/examples/full-ray-tracing/summary.rml new file mode 100644 index 00000000..78a511eb --- /dev/null +++ b/examples/full-ray-tracing/summary.rml @@ -0,0 +1,114 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/likelihood_BabyIAXO.rml b/examples/likelihood_BabyIAXO.rml deleted file mode 100644 index cf7fa5c2..00000000 --- a/examples/likelihood_BabyIAXO.rml +++ /dev/null @@ -1,28 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/likelihood_BabyIAXO_2.rml b/examples/likelihood_BabyIAXO_2.rml deleted file mode 100644 index 4a71524d..00000000 --- a/examples/likelihood_BabyIAXO_2.rml +++ /dev/null @@ -1,12 +0,0 @@ - - - - - - - - - - - - diff --git a/examples/metadata.rml b/examples/metadata.rml deleted file mode 100644 index 7e3c00c8..00000000 --- a/examples/metadata.rml +++ /dev/null @@ -1,48 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/axionGenerator.rml b/examples/old/axionGenerator.rml similarity index 100% rename from examples/axionGenerator.rml rename to examples/old/axionGenerator.rml diff --git a/examples/axionSpectrum.xml b/examples/old/axionSpectrum.xml similarity index 100% rename from examples/axionSpectrum.xml rename to examples/old/axionSpectrum.xml diff --git a/examples/axionlib.xml b/examples/old/axionlib.xml similarity index 100% rename from examples/axionlib.xml rename to examples/old/axionlib.xml diff --git a/examples/bField_BabyIAXO.rml b/examples/old/bField_BabyIAXO.rml similarity index 68% rename from examples/bField_BabyIAXO.rml rename to examples/old/bField_BabyIAXO.rml index 18ebac2e..eb06bc99 100644 --- a/examples/bField_BabyIAXO.rml +++ b/examples/old/bField_BabyIAXO.rml @@ -1,6 +1,5 @@ - - - - - + + + + diff --git a/examples/old/likelihood_BabyIAXO.rml b/examples/old/likelihood_BabyIAXO.rml new file mode 100644 index 00000000..10c0b7dd --- /dev/null +++ b/examples/old/likelihood_BabyIAXO.rml @@ -0,0 +1,24 @@ + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/old/likelihood_BabyIAXO_2.rml b/examples/old/likelihood_BabyIAXO_2.rml new file mode 100644 index 00000000..657e848d --- /dev/null +++ b/examples/old/likelihood_BabyIAXO_2.rml @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/examples/old/metadata.rml b/examples/old/metadata.rml new file mode 100644 index 00000000..280d5c08 --- /dev/null +++ b/examples/old/metadata.rml @@ -0,0 +1,39 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/opticsBench.rml b/examples/optics/opticsBench.rml similarity index 100% rename from examples/opticsBench.rml rename to examples/optics/opticsBench.rml diff --git a/examples/opticsBenchTrueW1.rml b/examples/optics/opticsBenchTrueW1.rml similarity index 80% rename from examples/opticsBenchTrueW1.rml rename to examples/optics/opticsBenchTrueW1.rml index d0969112..063b332f 100644 --- a/examples/opticsBenchTrueW1.rml +++ b/examples/optics/opticsBenchTrueW1.rml @@ -23,32 +23,25 @@ - + - - - - - + + + - - - - - + + + + + + + - - - - - - - + @@ -58,7 +51,7 @@ - + @@ -74,7 +67,6 @@ - @@ -92,20 +84,18 @@ - - - - - - - - + + + + + + diff --git a/examples/opticsPlots.rml b/examples/optics/opticsPlots.rml similarity index 100% rename from examples/opticsPlots.rml rename to examples/optics/opticsPlots.rml diff --git a/examples/opticsPlotsAllTransm.rml b/examples/optics/opticsPlotsAllTransm.rml similarity index 100% rename from examples/opticsPlotsAllTransm.rml rename to examples/optics/opticsPlotsAllTransm.rml diff --git a/examples/test.C b/examples/test.C deleted file mode 100644 index 0accd9a7..00000000 --- a/examples/test.C +++ /dev/null @@ -1,16 +0,0 @@ -void test() { - TFile f("Run_babyIAXO_1Volume_XY.root"); - // TFile f("Run_babyIAXO_1Volume_sphereIn.root"); - // TFile f("Run_babyIAXO_1Volume_sphereOut.root"); - TRestAnalysisTree* anaTree; - f.GetObject("AnalysisTree", anaTree); - // anaTree->Draw("axGen2_posY:axGen2_posZ"); - // anaTree->Scan("axGen2_posX"); - // anaTree->Scan("axGen2_probability"); - anaTree->Draw("axGen2_probability:axGen2_posY:axGen2_posX", "", "", 10000, - 0); // probability for the plan - // anaTree->Draw("axGen2_probability>0:axGen2_posY:axGen2_posX","","",10000,0); // To get the shape of the - // geometry - // anaTree->Draw("axGen2_probability:TMath::ATan(axGen2_posZ/axGen2_posY):TMath::ACos(axGen2_posX/15000)","","",10000,0); - // // probability for the sphere -} diff --git a/images/axionFieldPlots.png b/images/axionFieldPlots.png new file mode 100644 index 00000000..85eaa765 Binary files /dev/null and b/images/axionFieldPlots.png differ diff --git a/inc/TRestAxionEvent.h b/inc/TRestAxionEvent.h index 183f95da..56ff7c10 100644 --- a/inc/TRestAxionEvent.h +++ b/inc/TRestAxionEvent.h @@ -47,17 +47,8 @@ class TRestAxionEvent : public TRestEvent { /// Axion mass in eV Double_t fMass = 0.; //< - /// The effective magnetic field fixed by TRestAxionFieldPropagationProcess - Double_t fBSquared = 0; //< - - /// The effective conversion length fixed by TRestAxionFieldPropagationProcess - Double_t fLcoherence = 0; //< - /// It keeps track of efficiency introduced at different helioscope components - std::map fEfficiencies; - - /// We may use it to integrate a detector response inside each event - std::vector fResponse; //< + // std::map fEfficiencies; //< protected: public: @@ -76,10 +67,7 @@ class TRestAxionEvent : public TRestEvent { Double_t GetEnergy() { return fEnergy; } // returns value in keV Double_t GetMass() { return fMass * units("eV"); } // returns value in eV - Double_t GetBSquared() { return fBSquared; } - Double_t GetLConversion() { return fLcoherence; } - - Double_t GetEfficiency(std::string name) { return fEfficiencies[name]; } + // Double_t GetEfficiency(std::string name) { return fEfficiencies[name]; } void SetPosition(TVector3 pos) { fPosition = pos; } void SetPosition(Double_t x, Double_t y, Double_t z) { SetPosition(TVector3(x, y, z)); } @@ -98,10 +86,7 @@ class TRestAxionEvent : public TRestEvent { void SetEnergy(Double_t en) { fEnergy = en; } void SetMass(Double_t m) { fMass = m; } - void SetBSquared(Double_t b) { fBSquared = b; } - void SetLConversion(Double_t conv) { fLcoherence = conv; } - - void AddEfficiency(std::string name, Double_t value) { fEfficiencies[name] = value; } + // void AddEfficiency(std::string name, Double_t value) { fEfficiencies[name] = value; } virtual void Initialize(); @@ -112,6 +97,6 @@ class TRestAxionEvent : public TRestEvent { TRestAxionEvent(); ~TRestAxionEvent(); - ClassDef(TRestAxionEvent, 1); + ClassDef(TRestAxionEvent, 2); }; #endif diff --git a/inc/TRestAxionEventProcess.h b/inc/TRestAxionEventProcess.h index df6147c3..f91e615a 100644 --- a/inc/TRestAxionEventProcess.h +++ b/inc/TRestAxionEventProcess.h @@ -38,6 +38,9 @@ class TRestAxionEventProcess : public TRestEventProcess { /// The rotation angle with respect to X-axis Double_t fPitch = 0; + /// If enabled it will skip the end rotation that recovers the original axion trajectory direction + Bool_t fSkipEndProcessRotation = false; //! + protected: /// A pointer to the specific TRestAxionEvent TRestAxionEvent* fAxionEvent; //! @@ -47,6 +50,8 @@ class TRestAxionEventProcess : public TRestEventProcess { TVector3 GetCenter() const { return fCenter; } + void SkipEndProcessRotation(Bool_t value = true) { fSkipEndProcessRotation = value; } + public: RESTValue GetInputEvent() const override { return fAxionEvent; } RESTValue GetOutputEvent() const override { return fAxionEvent; } diff --git a/inc/TRestAxionField.h b/inc/TRestAxionField.h new file mode 100644 index 00000000..22662d63 --- /dev/null +++ b/inc/TRestAxionField.h @@ -0,0 +1,65 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +#ifndef _TRestAxionField +#define _TRestAxionField + +#include "TRestAxionBufferGas.h" + +//! A basic class to define analytical axion-photon conversion calculations for axion helioscopes +class TRestAxionField : public TObject { + private: + Bool_t fDebug = false; //! + + void Initialize(); + + /// A pointer to the buffer gas definition + TRestAxionBufferGas* fBufferGas = NULL; //! + + public: + Double_t BL(Double_t Bmag, Double_t Lcoh); + Double_t BLHalfSquared(Double_t Bmag, Double_t Lcoh); + + /// It enables/disables debug mode + void SetDebug(Bool_t v) { fDebug = v; } + + /// It assigns a gas buffer medium to the calculation + void AssignBufferGas(TRestAxionBufferGas* buffGas) { fBufferGas = buffGas; } + + /// It assigns a gas buffer medium to the calculation + void SetBufferGas(TRestAxionBufferGas* buffGas) { fBufferGas = buffGas; } + + Double_t GammaTransmissionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, + Double_t mg = 0, Double_t absLength = 0); + + Double_t GammaTransmissionProbability(std::vector Bmag, Double_t deltaL, Double_t Ea, + Double_t ma, Double_t mg = 0, Double_t absLength = 0); + + Double_t AxionAbsorptionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, + Double_t mg = 0, Double_t absLength = 0); + + TRestAxionField(); + ~TRestAxionField(); + + ClassDef(TRestAxionField, 2); +}; +#endif diff --git a/inc/TRestAxionFieldPropagationProcess.h b/inc/TRestAxionFieldPropagationProcess.h index b617705a..287645db 100644 --- a/inc/TRestAxionFieldPropagationProcess.h +++ b/inc/TRestAxionFieldPropagationProcess.h @@ -28,14 +28,27 @@ #include "TRestAxionEvent.h" #include "TRestAxionEventProcess.h" +#include "TRestAxionField.h" #include "TRestAxionMagneticField.h" #include "TRestPhysics.h" //! A process to introduce the magnetic field profile integration along the track class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess { private: + /// The differential length in mm used for the field integration + Double_t fIntegrationStep = 50; //< + + /// The additional length in mm that the converted photon propagates without magnetic field + Double_t fBufferGasAdditionalLength = 0; //< + /// A pointer to the magnetic field description stored in TRestRun - TRestAxionMagneticField* fField; //! + TRestAxionMagneticField* fMagneticField = nullptr; //! + + /// A pointer to TRestAxionField that implements probability calculations + TRestAxionField* fAxionField = nullptr; //! + + /// A pointer to TRestBufferGas given to TRestAxionField to perform calculations in a particular gas + TRestAxionBufferGas* fBufferGas = nullptr; //! void Initialize() override; @@ -47,6 +60,17 @@ class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess { TRestEvent* ProcessEvent(TRestEvent* eventInput) override; + /// It prints out the process parameters stored in the metadata structure + void PrintMetadata() override { + BeginPrintProcess(); + + RESTMetadata << "Integration step length : " << fIntegrationStep << " mm" << RESTendl; + RESTMetadata << "Buffer gas additional length : " << fBufferGasAdditionalLength * units("m") << " m" + << RESTendl; + + EndPrintProcess(); + } + /// Returns the name of this process const char* GetProcessName() const override { return "axionFieldPropagation"; } @@ -57,6 +81,6 @@ class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess { // Destructor ~TRestAxionFieldPropagationProcess(); - ClassDefOverride(TRestAxionFieldPropagationProcess, 1); + ClassDefOverride(TRestAxionFieldPropagationProcess, 2); }; #endif diff --git a/inc/TRestAxionLikelihood.h b/inc/TRestAxionLikelihood.h index afd8d231..5a0c70d6 100644 --- a/inc/TRestAxionLikelihood.h +++ b/inc/TRestAxionLikelihood.h @@ -26,7 +26,7 @@ #include #include "TRestAxionBufferGas.h" -#include "TRestAxionPhotonConversion.h" +#include "TRestAxionField.h" #include "TRestAxionSolarModel.h" #include "TRestAxionSpectrum.h" @@ -63,9 +63,9 @@ class TRestAxionLikelihood : public TRestMetadata { Double_t fLastStepDensity = 0.; //-> - TRestAxionPhotonConversion* fPhotonConversion; //! - TRestAxionBufferGas* fBufferGas; //! - TRestAxionSpectrum* fAxionSpectrum; //! + TRestAxionField* fAxionField; //! + TRestAxionBufferGas* fBufferGas; //! + TRestAxionSpectrum* fAxionSpectrum; //! /// Random number generator TRandom3* fRandom; //! diff --git a/inc/TRestAxionOpticsProcess.h b/inc/TRestAxionOpticsProcess.h index e79514d3..563025c4 100644 --- a/inc/TRestAxionOpticsProcess.h +++ b/inc/TRestAxionOpticsProcess.h @@ -57,9 +57,6 @@ class TRestAxionOpticsProcess : public TRestAxionEventProcess { TRestEvent* ProcessEvent(TRestEvent* evInput) override; - /// End of event process. Called directly after ProcessEvent() - void EndOfEventProcess(TRestEvent* evInput = nullptr) override; - void LoadConfig(std::string cfgFilename, std::string name = ""); /// It prints out the process parameters stored in the metadata structure diff --git a/inc/TRestAxionPhotonConversion.h b/inc/TRestAxionPhotonConversion.h deleted file mode 100644 index 5065272b..00000000 --- a/inc/TRestAxionPhotonConversion.h +++ /dev/null @@ -1,172 +0,0 @@ -/************************************************************************* - * 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. * - *************************************************************************/ - -#ifndef _TRestAxionPhotonConversion -#define _TRestAxionPhotonConversion - -#include "TRestAxionBufferGas.h" -#include "mpreal.h" - -/* -/// MOVED TO TRestAxionFieldPropagationProcess class -/// A structure to define the two components of a complex number using real precision. -/// To be used inside TRestAxionPhotonConversion. -struct ComplexReal { - /// The real part of the number - mpfr::mpreal real = 0; - - /// The imaginary part of the number - mpfr::mpreal img = 0; -}; -*/ -//! A basic class to define analytical axion-photon conversion calculations for axion helioscopes -class TRestAxionPhotonConversion : public TObject { - private: - /// A two component vector to store the complex EM field amplitude. - /// MOVED to TRestAxionFieldPropagationProcess - /// ComplexReal fAem; //! - - /// A two component vector to store the complex axion field amplitude. - /// MOVED to TRestAxionFieldPropagationProcess - /// ComplexReal faxion; //! - - Bool_t fDebug = false; //! - - void Initialize(); - - /// A pointer to the buffer gas definition - TRestAxionBufferGas* fBufferGas = NULL; //! - - /* - /// MOVED TO TRestAxionFieldPropagationProcess class - ///////////////////////////////////////////////////////////////////////// - // ----- Just a quick implementation of complex number operations ---- // - // ------------ including mpfr real precision arithmetics ------------ // - - /// A function to calculate complex number addition with real precision - ComplexReal ComplexAddition(const ComplexReal& a, const ComplexReal& b) { - ComplexReal c; - - c.real = a.real + b.real; - c.img = a.img + b.img; - - return c; - } - - /// A function to calculate complex number substraction with real precision - ComplexReal ComplexSubstraction(const ComplexReal& a, const ComplexReal& b) { - ComplexReal c; - - c.real = a.real - b.real; - c.img = a.img - b.img; - - return c; - } - - /// A function to calculate complex number product with real precision - ComplexReal ComplexProduct(const ComplexReal& a, const ComplexReal& b) { - ComplexReal c; - - c.real = a.real * b.real - a.img * b.img; - c.img = a.real * b.img + a.img * b.real; - - return c; - } - - /// A function to calculate complex number product by a value with real precision - ComplexReal ComplexProduct(const mpfr::mpreal& value, const ComplexReal& a) { - ComplexReal c; - - c.real = value * a.real; - c.img = value * a.img; - - return c; - } - - /// A function to calculate complex number cocient with real precision - ComplexReal ComplexCocient(const ComplexReal& a, const ComplexReal& b) { - ComplexReal c = ComplexConjugate(b); - c = ComplexProduct(a, c); - - mpfr::mpreal norm = 1. / Norm2(b); - - c = ComplexProduct(norm, c); - - return c; - } - - /// A function to calculate complex conjugate with real precision - ComplexReal ComplexConjugate(const ComplexReal& a) { - ComplexReal c; - - c.real = a.real; - c.img = -a.img; - - return c; - } - - /// A function to calculate the norm squared from a complex number with real precision - mpfr::mpreal Norm2(const ComplexReal& a) { - mpfr::mpreal result = a.real * a.real + a.img * a.img; - return result; - } - - /// A function to calculate complex number product by a value with real precision - ComplexReal SetComplexReal(const mpfr::mpreal& r, const mpfr::mpreal& i) { - ComplexReal c; - - c.real = r; - c.img = i; - - return c; - } - */ - public: - Double_t BL(Double_t Bmag, Double_t Lcoh); - Double_t BLHalfSquared(Double_t Bmag, Double_t Lcoh); - - /// It enables/disables debug mode - void SetDebug(Bool_t v) { fDebug = v; } - - /// It assigns a gas buffer medium to the calculation - void AssignBufferGas(TRestAxionBufferGas* buffGas) { fBufferGas = buffGas; } - - /// It assigns a gas buffer medium to the calculation - void SetBufferGas(TRestAxionBufferGas* buffGas) { fBufferGas = buffGas; } - - Double_t GammaTransmissionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, - Double_t mg = 0, Double_t absLength = 0); - - Double_t AxionAbsorptionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, - Double_t mg = 0, Double_t absLength = 0); - - /// Commented because it uses ComplexReal structure that is moved to TRestAxionFieldPropagationProcess - /// class - /// void PropagateAxion(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, Double_t mg = 0, - /// Double_t absLength = 0); - - TRestAxionPhotonConversion(); - ~TRestAxionPhotonConversion(); - - ClassDef(TRestAxionPhotonConversion, 2); -}; -#endif diff --git a/inc/TRestAxionTrueWolterOptics.h b/inc/TRestAxionTrueWolterOptics.h index 10f55073..505d6465 100644 --- a/inc/TRestAxionTrueWolterOptics.h +++ b/inc/TRestAxionTrueWolterOptics.h @@ -128,7 +128,9 @@ class TRestAxionTrueWolterOptics : public TRestAxionOptics { } void SetMirror() override { - fCurrentMirror = fEntranceRingsMask->GetRegion(fEntrancePosition.X(), fEntrancePosition.Y()); + Double_t x = fEntrancePosition.X(); + Double_t y = fEntrancePosition.Y(); + fCurrentMirror = fEntranceRingsMask->GetRegion(x, y); } TPad* DrawMirrors() override; diff --git a/inc/TRestAxionWolterOptics.h b/inc/TRestAxionWolterOptics.h index dbf3e119..83be3ccb 100644 --- a/inc/TRestAxionWolterOptics.h +++ b/inc/TRestAxionWolterOptics.h @@ -128,7 +128,9 @@ class TRestAxionWolterOptics : public TRestAxionOptics { } void SetMirror() override { - fCurrentMirror = fEntranceRingsMask->GetRegion(fEntrancePosition.X(), fEntrancePosition.Y()); + Double_t x = fEntrancePosition.X(); + Double_t y = fEntrancePosition.Y(); + fCurrentMirror = fEntranceRingsMask->GetRegion(x, y); } TPad* DrawMirrors() override; diff --git a/pipeline/physics/AxionPhysicsValidation.C b/pipeline/physics/AxionPhysicsValidation.C new file mode 100644 index 00000000..78a1be6f --- /dev/null +++ b/pipeline/physics/AxionPhysicsValidation.C @@ -0,0 +1,180 @@ +#include +#include + +Int_t AxionPhysicsValidation() { + std::cout << endl; + + std::cout << "Evaluating magnetic field for BabyIAXO" << std::endl; + std::cout << "--------------------------------------" << std::endl; + + TRestAxionMagneticField magneticField("fields.rml", "babyIAXO"); + + std::cout << "Field at (0,0,0) : " << magneticField.GetMagneticField(0, 0, 0).Y() << std::endl; + std::cout << "Field at (0,0,-4000) : " << magneticField.GetMagneticField(0, 0, -4000).Y() << std::endl; + std::cout << "Field at (0,300,4000) : " << magneticField.GetMagneticField(0, 300, 4000).Y() << std::endl; + + if (std::round(magneticField.GetMagneticField(0, 0, 0).Y() * 1000) != -2007) { + std::cout << "Field at (0,0,0) is not -2.07 T!" << std::endl; + return 10; + } + + if (std::round(magneticField.GetMagneticField(0, 0, -4000).Y() * 1000) != -1700) { + std::cout << "Field at (0,0,-4000) is not -1.7 T!" << std::endl; + return 11; + } + + if (std::round(magneticField.GetMagneticField(0, 300, 4000).Y() * 1000) != -1290) { + std::cout << "Field at (0,0,4000) is not -1.29 T!" << std::endl; + return 12; + } + + std::cout << endl; + std::cout << "--> Magnetic field evaluation succeeded!" << std::endl; + std::cout << endl; + + std::cout << "Evaluating dummy gas mixtures" << std::endl; + std::cout << "-----------------------------" << std::endl; + + TRestAxionBufferGas helium("bufferGases.rml", "helium"); + + std::cout << "Gas name : " << helium.GetName() << std::endl; + std::cout << "---------------" << std::endl; + std::cout << "Number of gases : " << helium.GetNumberOfGases() << std::endl; + std::cout << "Photon mass at 3keV : " << helium.GetPhotonMass(3) << " eV" << std::endl; + std::cout << "Absorption lenght (cm-1) at 3keV : " << helium.GetPhotonAbsorptionLength(3) << std::endl; + + if (std::round(1.e6 * helium.GetPhotonMass(3)) != 1017) { + std::cout << "Photon mass is not 1.017 meV!" << std::endl; + return 21; + } + + if (std::round(1.e12 * helium.GetPhotonAbsorptionLength(3)) != 4720) { + std::cout << "Absorption length is not 4.72e-9 cm-1!" << std::endl; + return 22; + } + + std::cout << std::endl; + + TRestAxionBufferGas mixture("bufferGases.rml", "HeNeMixture"); + + std::cout << "Gas name : " << mixture.GetName() << std::endl; + std::cout << "---------------" << std::endl; + std::cout << "Number of gases : " << mixture.GetNumberOfGases() << std::endl; + std::cout << "Photon mass at 3keV : " << mixture.GetPhotonMass(3) << " eV" << std::endl; + std::cout << "Absorption lenght (cm-1) at 3keV : " << mixture.GetPhotonAbsorptionLength(3) << std::endl; + + if (std::round(1.e3 * mixture.GetPhotonMass(3)) != 651) { + std::cout << "Photon mass is not 0.651 eV!" << std::endl; + return 31; + } + + if (std::round(1.e3 * mixture.GetPhotonAbsorptionLength(3)) != 397) { + std::cout << "Absorption length is not 0.397 cm-1!" << std::endl; + return 32; + } + + std::cout << endl; + std::cout << "--> Buffer gas evaluation succeeded!" << std::endl; + std::cout << endl; + + //// Retrieving the magnetic field profile for a dummy trajectory. + //// We force it to be strongly un-aligned with z + + TVector3 position1(-300, 600, -6000); + TVector3 position2(600, -600, 6000); + + TVector3 direction = position2 - position1; + direction.Unit().Print(); + + std::vector v = magneticField.GetFieldBoundaries(position1, direction.Unit()); + if (v.size() != 2) { + std::cout << "The track does not traverse the magnetic volume" << std::endl; + return 100; + } + + std::cout << "In position. X: " << v[0].X() << " Y: " << v[0].Y() << " Z: " << v[0].Z() << std::endl; + std::cout << "Out position. X: " << v[1].X() << " Y: " << v[1].Y() << " Z: " << v[1].Z() << std::endl; + + TVector3 in = v[0]; + TVector3 out = v[1]; + + if (std::round(v[0].X()) != -101 || std::round(v[0].Y()) != 335 || std::round(v[0].Z()) != -3350) { + std::cout << "Wrong boundary determination. The input position is wrong!" << std::endl; + return 101; + } + + if (std::round(v[1].X()) != 293 || std::round(v[1].Y()) != -191 || std::round(v[1].Z()) != 1910) { + std::cout << "Wrong boundary determination. The output position is wrong!" << std::endl; + return 102; + } + + Double_t dl = 50; // 5cm + std::vector bProfile = magneticField.GetTransversalComponentAlongPath(v[0], v[1], dl); + + // for (auto& x : bProfile) std::cout << x << " "; + // std::cout << std::endl; + + std::cout << "Field elements (5cm step) : " << bProfile.size() << std::endl; + if (bProfile.size() != 107) { + std::cout << "Number of field elements is not the right value!" << std::endl; + return 103; + } + + std::cout << "Coherence length : " << (bProfile.size() - 1) * 50 << " mm" << std::endl; + + Double_t fieldAverage = std::accumulate(bProfile.begin(), bProfile.end(), 0.0) / bProfile.size(); + + std::cout << "Field average (T): " << fieldAverage << std::endl; + + if (std::round(fieldAverage * 1000.) != 1901) { + std::cout << "Field average is not 1.901 T!" << std::endl; + return 104; + } + + std::cout << endl; + std::cout << "--> Track boundaries and field along path evaluation succeeded!" << std::endl; + std::cout << endl; + + TRestAxionField axionField; + axionField.AssignBufferGas(&helium); + + std::cout << "BL(8.8T, 9.26m) : " << axionField.BL(8.8, 9260) << std::endl; + std::cout << "0.5 x (BL)^2 : " << axionField.BLHalfSquared(8.8, 9260) << std::endl; + + if (std::round(axionField.BL(8.8, 9260) * 1e10) != 81) { + std::cout << "BL calculation is not 8.1e-9!" << std::endl; + return 201; + } + + if (std::round(axionField.BLHalfSquared(8.8, 9260) * 1e20) != 1627) { + std::cout << "(BL)^2/2 calculation is not 1.627e-17!" << std::endl; + return 202; + } + + Double_t B = 8.8; // T + Double_t L = 9260.; // mm + Double_t ma = 1.e-2; // eV + Double_t Ea = 3.0; // keV + + std::cout << "Probability (CAST. Helium. ma:10meV - Ea:3keV): " + << axionField.GammaTransmissionProbability(B, L, Ea, ma) << std::endl; + + if (std::round(axionField.GammaTransmissionProbability(B, L, Ea, ma) * 1.e20) != 1547) { + std::cout << "CAST Probability. Wrong axion-photon conversion probability" << std::endl; + return 301; + } + + Double_t prob = axionField.GammaTransmissionProbability(bProfile, 50, Ea, ma); + std::cout << "Probability BabyIAXO : " << prob << std::endl; + + if (std::round(prob * 1.e22) != 2482) { + std::cout << "BabyIAXO Probability. Wrong axion-photon conversion probability" << std::endl; + return 401; + } + + std::cout << endl; + std::cout << "--> Axion-photon probability calculations succeeded!" << std::endl; + std::cout << endl; + + return 0; +} diff --git a/pipeline/ray-tracing/axion-field/Validate.C b/pipeline/ray-tracing/axion-field/Validate.C new file mode 100644 index 00000000..598eca5c --- /dev/null +++ b/pipeline/ray-tracing/axion-field/Validate.C @@ -0,0 +1,51 @@ +#include +#include + +Double_t probTolerance = 1.e-21; +Double_t fieldTolerance = 0.01; + +Int_t Validate(Double_t prob = 5.9978e-19, Double_t fieldAverage = 1.631) { + TRestRun* run = new TRestRun("AxionPhotonProbability.root"); + + if (run->GetEntries() != 100) { + std::cout << "Error. Number of entries is not 10000!" << std::endl; + + return 1; + } + + run->GetAnalysisTree()->Draw("axionPhoton_probability", "axionPhoton_probability"); + + TH1D* h = (TH1D*)run->GetAnalysisTree()->GetHistogram(); + if (h == nullptr) { + std::cout << "Problems generating probability histogram" << std::endl; + return 2; + } + + Double_t probability = h->Integral() / run->GetEntries(); + std::cout << "Average axion-photon probability: " << probability << std::endl; + + if (probability < prob - probTolerance || probability > prob + probTolerance) { + std::cout << "Wrong axion-photon probability!" << std::endl; + return 3; + } + + run->GetAnalysisTree()->Draw("axionPhoton_fieldAverage", "axionPhoton_fieldAverage"); + + TH1D* h2 = (TH1D*)run->GetAnalysisTree()->GetHistogram(); + if (h2 == nullptr) { + std::cout << "Problems generating field average histogram" << std::endl; + return 4; + } + + Double_t field = h2->Integral() / run->GetEntries(); + std::cout << "Average magnetic field: " << field << " T" << std::endl; + + if (field < fieldAverage - fieldTolerance || field > fieldAverage + fieldTolerance) { + std::cout << "Wrong average field!" << std::endl; + return 5; + } + + delete run; + + return 0; +} diff --git a/pipeline/ray-tracing/axion-field/photonConversion.rml b/pipeline/ray-tracing/axion-field/photonConversion.rml new file mode 100644 index 00000000..1d176dc4 --- /dev/null +++ b/pipeline/ray-tracing/axion-field/photonConversion.rml @@ -0,0 +1,33 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/pipeline/ray-tracing/axion-field/plots.rml b/pipeline/ray-tracing/axion-field/plots.rml new file mode 100644 index 00000000..bdc704ff --- /dev/null +++ b/pipeline/ray-tracing/axion-field/plots.rml @@ -0,0 +1,50 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/pipeline/ray-tracing/ValidateXMM.C b/pipeline/ray-tracing/optics/ValidateXMM.C similarity index 92% rename from pipeline/ray-tracing/ValidateXMM.C rename to pipeline/ray-tracing/optics/ValidateXMM.C index 2cca73f4..8bf1ebcf 100644 --- a/pipeline/ray-tracing/ValidateXMM.C +++ b/pipeline/ray-tracing/optics/ValidateXMM.C @@ -13,6 +13,9 @@ Int_t ValidateXMM(std::string fname) { run->GetAnalysisTree()->Draw("optics_efficiency", "optics_efficiency"); + run->GetEntry(10); + run->GetAnalysisTree()->PrintObservables(); + TH1D* h = (TH1D*)run->GetAnalysisTree()->GetHistogram(); if (h == nullptr) { diff --git a/pipeline/validateMacros.py b/pipeline/validateMacros.py new file mode 100755 index 00000000..005d339b --- /dev/null +++ b/pipeline/validateMacros.py @@ -0,0 +1,21 @@ +import os, sys + +os.system("restRoot -b -q -m 1 > error.log 2>&1") + +with open('error.log') as f: + lines = f.readlines() + +for line in lines: + if line.find("warning") >= 0 or line.find("Warning") >= 0: + print("Find warnings when launching restRoot!") + print(lines) + sys.exit(1) + if line.find("error") >= 0 or line.find("Error") >= 0: + print("Find error messages when launching restRoot!") + print(lines) + sys.exit(2) + +print("No warning found when loading macros") +print("Succeed") + +sys.exit(0) diff --git a/src/TRestAxionAnalysisProcess.cxx b/src/TRestAxionAnalysisProcess.cxx index 4a53e24d..df968b5b 100644 --- a/src/TRestAxionAnalysisProcess.cxx +++ b/src/TRestAxionAnalysisProcess.cxx @@ -21,7 +21,17 @@ *************************************************************************/ ////////////////////////////////////////////////////////////////////////// -/// TRestAxionAnalysisProcess TOBE documented +/// TRestAxionAnalysisProcess adds to the analysis tree information about +/// the particle state. +/// +/// Observable values defined by this process: +/// - **energy**: The energy in keV. +/// - **posX**: The X-position in mm. +/// - **posY**: The Y-position in mm. +/// - **posZ**: The Z-position in mm. +/// - **thetaAngle**: The angle in radians respect to the propagation axis, defined in Z. +/// - **phiAngle**: The angle moving along the XY-plane. +/// - **R**: The distance to the Z-axis in mm. /// /// The axion is generated with intensity proportional to g_ag = 1.0 x g10 ///-------------------------------------------------------------------------- @@ -114,12 +124,17 @@ TRestEvent* TRestAxionAnalysisProcess::ProcessEvent(TRestEvent* evInput) { SetObservableValue("energy", fAxionEvent->GetEnergy()); - SetObservableValue("posX", fAxionEvent->GetPosition().X()); - SetObservableValue("posY", fAxionEvent->GetPosition().Y()); + Double_t x = fAxionEvent->GetPosition().X(); + Double_t y = fAxionEvent->GetPosition().Y(); + SetObservableValue("posX", x); + SetObservableValue("posY", y); SetObservableValue("posZ", fAxionEvent->GetPosition().Z()); - // SetObservableValue("B2", fAxionEvent->GetBSquared()); - // SetObservableValue("Lcoh", fAxionEvent->GetLConversion()); + Double_t r = TMath::Sqrt(x * x + y * y); + SetObservableValue("R", r); + + SetObservableValue("thetaAngle", fAxionEvent->GetDirection().Theta()); + SetObservableValue("phiAngle", fAxionEvent->GetDirection().Phi()); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fAxionEvent->PrintEvent(); diff --git a/src/TRestAxionBufferGas.cxx b/src/TRestAxionBufferGas.cxx index 6ef19440..c19a7958 100644 --- a/src/TRestAxionBufferGas.cxx +++ b/src/TRestAxionBufferGas.cxx @@ -142,8 +142,13 @@ void TRestAxionBufferGas::InitFromConfigFile() { auto gasDefinition = GetElement("gas"); while (gasDefinition) { TString gasName = GetFieldValue("name", gasDefinition); - Double_t gasDensity = GetDblParameterWithUnits("density", gasDefinition) * units("g/cm3"); - SetGasDensity(gasName, gasDensity); + if (gasName.Contains("+")) { + TString gasDensity = GetFieldValue("density", gasDefinition); + SetGasMixture(gasName, gasDensity); + } else { + Double_t gasDensity = GetDblParameterWithUnits("density", gasDefinition); + SetGasDensity(gasName, gasDensity); + } gasDefinition = GetNextElement(gasDefinition); } @@ -153,6 +158,8 @@ void TRestAxionBufferGas::InitFromConfigFile() { /////////////////////////////////////////////// /// \brief It adds a new gas component to the mixture. If it already exists it will update its density. /// +/// Density must be given in standard REST units: kg/mm^3 +/// void TRestAxionBufferGas::SetGasDensity(TString gasName, Double_t density) { Int_t gasIndex = FindGasIndex(gasName); @@ -176,7 +183,7 @@ void TRestAxionBufferGas::SetGasMixture(TString gasMixture, TString gasDensities if (names.size() == densities.size()) { for (int n = 0; n < names.size(); n++) { Double_t density = GetValueInRESTUnits(densities[n]); - SetGasDensity(names[n], density * units("g/cm3")); + SetGasDensity(names[n], density); } } else { this->SetError("SetGasMixture. Number of gases does not match the densities!"); @@ -340,7 +347,8 @@ Double_t TRestAxionBufferGas::GetFormFactor(TString gasName, Double_t energy) { Double_t TRestAxionBufferGas::GetPhotonAbsorptionLength(Double_t energy) { Double_t attLength = 0; for (unsigned int n = 0; n < fBufferGasName.size(); n++) - attLength += fBufferGasDensity[n] * GetAbsorptionCoefficient(fBufferGasName[n], energy); + attLength += + fBufferGasDensity[n] * units("g/cm^3") * GetAbsorptionCoefficient(fBufferGasName[n], energy); return attLength; } @@ -381,7 +389,8 @@ Double_t TRestAxionBufferGas::GetPhotonMass(double en) { RESTError << "W value must be defined in TRestAxionBufferGas::GetPhotonMass" << RESTendl; RESTError << "This gas will not contribute to the calculation of the photon mass!" << RESTendl; } else { - photonMass += fBufferGasDensity[n] * GetFormFactor(fBufferGasName[n], en) / W_value; + photonMass += + fBufferGasDensity[n] * units("g/cm^3") * GetFormFactor(fBufferGasName[n], en) / W_value; } } @@ -497,11 +506,14 @@ void TRestAxionBufferGas::PrintMetadata() { if (fBufferGasName.size() == 0) { RESTMetadata << "Buffer medium is vacuum" << RESTendl; } else { - RESTMetadata << "Buffer gases defined : " << RESTendl; + RESTMetadata << "Photon mass at 4keV : " << this->GetPhotonMass(4.) << " eV" << RESTendl; + RESTMetadata << " " << RESTendl; + RESTMetadata << "Buffer gases inside mixture : " << RESTendl; RESTMetadata << "---------------------------" << RESTendl; for (unsigned int n = 0; n < fBufferGasName.size(); n++) { RESTMetadata << " Gas name : " << fBufferGasName[n] << RESTendl; - RESTMetadata << " Gas density : " << fBufferGasDensity[n] << " g/cm3" << RESTendl; + RESTMetadata << " Gas density : " << fBufferGasDensity[n] * units("g/cm^3") << " g/cm3" + << RESTendl; RESTMetadata << " Form factor energy range : ( " << fFactorEnergy[n][0] << ", " << fFactorEnergy[n].back() << " ) keV" << RESTendl; RESTMetadata << " Absorption energy range : ( " << fAbsEnergy[n][0] << ", " diff --git a/src/TRestAxionEvent.cxx b/src/TRestAxionEvent.cxx index 71ef12e8..ec262244 100644 --- a/src/TRestAxionEvent.cxx +++ b/src/TRestAxionEvent.cxx @@ -163,6 +163,5 @@ void TRestAxionEvent::PrintEvent() { << endl; cout << "Direction : ( " << fDirection.X() << ", " << fDirection.Y() << ", " << fDirection.Z() << " )" << endl; - cout << "B^2 : " << fBSquared << " T^2" << endl; cout << endl; } diff --git a/src/TRestAxionEventProcess.cxx b/src/TRestAxionEventProcess.cxx index b303a12a..12987f9e 100644 --- a/src/TRestAxionEventProcess.cxx +++ b/src/TRestAxionEventProcess.cxx @@ -98,7 +98,7 @@ void TRestAxionEventProcess::EndOfEventProcess(TRestEvent* evInput) { << RESTendl; RESTDebug << " ---- " << RESTendl; fAxionEvent->Translate(TVector3(fCenter.X(), fCenter.Y(), fCenter.Z())); - fAxionEvent->RotateXY(fCenter, fPitch, fYaw); + if (!fSkipEndProcessRotation) fAxionEvent->RotateXY(fCenter, fPitch, fYaw); RESTDebug << "EoEP: Final Position. X: " << fAxionEvent->GetPosition().X() << " Y: " << fAxionEvent->GetPosition().Y() << " Z: " << fAxionEvent->GetPosition().Z() << RESTendl; @@ -117,7 +117,7 @@ void TRestAxionEventProcess::EndOfEventProcess(TRestEvent* evInput) { void TRestAxionEventProcess::BeginPrintProcess() { TRestEventProcess::BeginPrintProcess(); - RESTMetadata << "Center: (" << fCenter.X() << ", " << fCenter.Y() << ", " << fCenter.Z() << ")" + RESTMetadata << "Center: (" << fCenter.X() << ", " << fCenter.Y() << ", " << fCenter.Z() << ") mm" << RESTendl; RESTMetadata << "Yaw angle (Y-axis): " << fYaw * units("degrees") << " degrees" << RESTendl; RESTMetadata << "Pitch angle (X-axis): " << fPitch * units("degrees") << " degrees" << RESTendl; diff --git a/src/TRestAxionPhotonConversion.cxx b/src/TRestAxionField.cxx similarity index 60% rename from src/TRestAxionPhotonConversion.cxx rename to src/TRestAxionField.cxx index b00e97d9..1514b285 100644 --- a/src/TRestAxionPhotonConversion.cxx +++ b/src/TRestAxionField.cxx @@ -21,12 +21,12 @@ *************************************************************************/ ////////////////////////////////////////////////////////////////////////// -/// TRestAxionPhotonConversion is a class used to calculate the axion-photon mixing +/// TRestAxionField is a class used to calculate the axion-photon mixing /// and determine the probability of the particle being in an axion or photon /// state. /// /// A peculiarity from this class is that it encapsulates internally the high -/// precision calculations using the real precisions types from library mpreal. +/// precision calculations using the real precisions types using TRestComplex. /// It is known that double precision is not good enough in some scenarios. /// ///-------------------------------------------------------------------------- @@ -35,45 +35,48 @@ /// /// History of developments: /// -/// 2019-March: First concept and implementation of TRestAxionPhotonConversion class. +/// 2019-March: First concept and implementation of TRestAxionField class. /// Javier Galan /// -/// \class TRestAxionPhotonConversion +/// \class TRestAxionField /// \author Javier Galan /// ///
/// -#include "TRestAxionPhotonConversion.h" +#include "TRestAxionField.h" #include -#include "TComplex.h" #include "TH1F.h" -using namespace std; -using namespace REST_Physics; +#ifdef USE_MPFR +#include "TRestComplex.h" +#endif + +#include -// Better we keep specifying mpfr:: explicitily -// using mpfr::mpreal; +using namespace std; -ClassImp(TRestAxionPhotonConversion); +ClassImp(TRestAxionField); /////////////////////////////////////////////// /// \brief Default constructor /// -TRestAxionPhotonConversion::TRestAxionPhotonConversion() { Initialize(); } +TRestAxionField::TRestAxionField() { Initialize(); } /////////////////////////////////////////////// /// \brief Default destructor /// -TRestAxionPhotonConversion::~TRestAxionPhotonConversion() {} +TRestAxionField::~TRestAxionField() {} /////////////////////////////////////////////// -/// \brief Initialization of TRestAxionPhotonConversion class +/// \brief Initialization of TRestAxionField class /// /// It sets the default real precision to be used with mpfr types. Now it is 30 digits. /// So that we can still calculate numbers such as : 1.0 - 1.e-30 /// -void TRestAxionPhotonConversion::Initialize() { - mpfr::mpreal::set_default_prec(mpfr::digits2bits(30)); +void TRestAxionField::Initialize() { +#ifdef USE_MPFR + TRestComplex::SetPrecision(30); +#endif fBufferGas = NULL; @@ -88,10 +91,10 @@ void TRestAxionPhotonConversion::Initialize() { /// `Lcoh` should be expressed in `mm`, and `Bmag` in `T`. /// The result will be given for an axion-photon coupling of 10^{-10} GeV^{-1} /// -double TRestAxionPhotonConversion::BL(Double_t Bmag, Double_t Lcoh) { +double TRestAxionField::BL(Double_t Bmag, Double_t Lcoh) { Double_t lengthInMeters = Lcoh / 1000.; - Double_t tm = lightSpeed / naturalElectron * 1.0e-9; // GeV + Double_t tm = REST_Physics::lightSpeed / REST_Physics::naturalElectron * 1.0e-9; // GeV Double_t sol = lengthInMeters * Bmag * tm; sol = sol * 1.0e-10; @@ -104,11 +107,11 @@ double TRestAxionPhotonConversion::BL(Double_t Bmag, Double_t Lcoh) { /// `Lcoh` should be expressed in `mm`, and `Bmag` in `T`. /// The result will be given for an axion-photon coupling of 10^{-10} GeV^{-1} /// -double TRestAxionPhotonConversion::BLHalfSquared(Double_t Bmag, Double_t Lcoh) // (BL/2)**2 +double TRestAxionField::BLHalfSquared(Double_t Bmag, Double_t Lcoh) // (BL/2)**2 { Double_t lengthInMeters = Lcoh / 1000.; - Double_t tm = lightSpeed / naturalElectron * 1.0e-9; // gev + Double_t tm = REST_Physics::lightSpeed / REST_Physics::naturalElectron * 1.0e-9; // gev Double_t sol = lengthInMeters * Bmag * tm / 2; sol = sol * sol * 1.0e-20; @@ -119,8 +122,9 @@ double TRestAxionPhotonConversion::BLHalfSquared(Double_t Bmag, Double_t Lcoh) /// \brief Performs the calculation of axion-photon conversion probability using directly /// equation (11) from van Bibber, Phys Rev D Part Fields. 1989. /// -/// m_gamma will be obtainned from buffer gas definition. If no buffer gas has been assigned the medium -/// will be assumed to be vacuum. +/// If m_gamma (mg) is not given as an argument, i.e. it is equal to zero, then m_gamma +/// will be obtainned from the buffer gas definition. If no buffer gas has been assigned +/// then the medium will be assumed to be vacuum. /// /// Ea in keV, ma in eV, mgamma in eV, Lcoh in mm, Bmag in T /// @@ -128,9 +132,15 @@ double TRestAxionPhotonConversion::BLHalfSquared(Double_t Bmag, Double_t Lcoh) /// /// The returned value is given for g_ag = 10^-10 GeV-1 /// -Double_t TRestAxionPhotonConversion::GammaTransmissionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, - Double_t ma, Double_t mg, - Double_t absLength) { +Double_t TRestAxionField::GammaTransmissionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, + Double_t mg, Double_t absLength) { +#ifndef USE_MPFR + RESTWarning + << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFR=ON to your REST compilation" + << RESTendl; + RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; + return 0; +#else mpfr::mpreal axionMass = ma; mpfr::mpreal cohLength = Lcoh / 1000.; // Default REST units are mm; @@ -139,7 +149,7 @@ Double_t TRestAxionPhotonConversion::GammaTransmissionProbability(Double_t Bmag, if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea); RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " TRestAxionPhotonConversion::GammaTransmissionProbability. Parameter summary" << RESTendl; + RESTDebug << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl; RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; RESTDebug << " Axion mass : " << ma << " eV" << RESTendl; RESTDebug << " Axion energy : " << Ea << " keV" << RESTendl; @@ -150,7 +160,7 @@ Double_t TRestAxionPhotonConversion::GammaTransmissionProbability(Double_t Bmag, if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(Bmag, Lcoh); mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / Ea / 1000.0; - mpfr::mpreal l = cohLength * PhMeterIneV; + mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; mpfr::mpreal phi = q * l; mpfr::mpreal Gamma = absLength; @@ -184,6 +194,125 @@ Double_t TRestAxionPhotonConversion::GammaTransmissionProbability(Double_t Bmag, RESTDebug << "Axion-photon transmission probability : " << sol << RESTendl; return sol; +#endif +} + +/////////////////////////////////////////////// +/// \brief Performs the calculation of axion-photon conversion probability using directly +/// equation (28) from J. Redondo and A. Ringwald, Light shinning through walls. +/// https://arxiv.org/pdf/1011.3741.pdf +/// +/// If m_gamma (mg) is not given as an argument, i.e. it is equal to zero, then m_gamma +/// will be obtainned from the buffer gas definition. If no buffer gas has been assigned +/// then the medium will be assumed to be vacuum. +/// +/// Ea in keV, ma in eV, mgamma in eV, deltaL in mm, Bmag in T +/// +/// mg in eV, absLength in cm-1 +/// +/// The returned value is given for g_ag = 10^-10 GeV-1 +/// +/// \note The density is for the moment homogeneous. We would need to implemnent a double integral +/// to solve the problem with a density profile. TOBE implemented in a new method if needed, where +/// Gamma is not constant and \integral{q(z)} is integrated at each step. +/// +Double_t TRestAxionField::GammaTransmissionProbability(std::vector Bmag, Double_t deltaL, + Double_t Ea, Double_t ma, Double_t mg, + Double_t absLength) { +#ifndef USE_MPFR + RESTWarning + << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFR=ON to your REST compilation" + << RESTendl; + RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; + return 0; +#else + mpfr::mpreal axionMass = ma; + + // Default REST units are mm. We express cohLength in m. + Double_t Lcoh = (Bmag.size() - 1) * deltaL; // in mm + Double_t cohLength = Lcoh / 1000.; // in m + + mpfr::mpreal photonMass = mg; + + if (mg == 0 && fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea); + + Double_t fieldAverage = 0; + if (Bmag.size() > 0) fieldAverage = std::accumulate(Bmag.begin(), Bmag.end(), 0.0) / Bmag.size(); + + RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; + RESTDebug << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl; + RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; + RESTDebug << " Axion mass : " << ma << " eV" << RESTendl; + RESTDebug << " Axion energy : " << Ea << " keV" << RESTendl; + RESTDebug << " Lcoh : " << cohLength << " mm" << RESTendl; + RESTDebug << " Bmag average : " << fieldAverage << " T" << RESTendl; + RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; + + // In vacuum + if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(fieldAverage, Lcoh); + + mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / Ea / 1000.0; + mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; + mpfr::mpreal phi = q * l; + + mpfr::mpreal Gamma = absLength; + if (absLength == 0 && fBufferGas) Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 + mpfr::mpreal GammaL = Gamma * cohLength * 100; + + if (fDebug) { + RESTDebug << "+------------------------+" << RESTendl; + RESTDebug << " Intermediate calculations" << RESTendl; + RESTDebug << " q : " << q << " eV" << RESTendl; + RESTDebug << " l : " << l << " eV-1" << RESTendl; + RESTDebug << " phi : " << phi << RESTendl; + RESTDebug << "Gamma : " << Gamma << RESTendl; + RESTDebug << "GammaL : " << GammaL << RESTendl; + RESTDebug << "+------------------------+" << RESTendl; + } + + mpfr::mpreal MFactor = phi * phi + GammaL * GammaL / 4.0; + MFactor = 1.0 / MFactor; + + if (fDebug) { + RESTDebug << "Mfactor : " << MFactor << RESTendl; + RESTDebug << "(BL/2)^2 : " << BLHalfSquared(fieldAverage, Lcoh) << RESTendl; + RESTDebug << "cos(phi) : " << cos(phi) << RESTendl; + RESTDebug << "Exp(-GammaL) : " << exp(-GammaL) << RESTendl; + } + + Double_t deltaIneV = deltaL / 1000. * REST_Physics::PhMeterIneV; + + /// We integrate following the Midpoint rule method. (Other potential options : Trapezoidal, Simpsons) + TRestComplex sum(0, 0); + for (unsigned int n = 0; n < Bmag.size() - 1; n++) { + Double_t Bmiddle = 0.5 * (Bmag[n] + Bmag[n + 1]); + + Double_t lStepIneV = ((double)n + 0.5) * deltaIneV; + Double_t lStepInCm = ((double)n + 0.5) * deltaL / 10.; + + TRestComplex qC(0, -q * lStepIneV); + qC = TRestComplex::Exp(qC); + + TRestComplex gC(0.5 * Gamma * lStepInCm, 0); + gC = TRestComplex::Exp(gC); + + TRestComplex integrand = Bmiddle * deltaL * gC * qC; // The integrand is in T by mm + + sum += integrand; + } + + mpfr::mpreal sol = exp(-GammaL) * sum.Rho2() * BLHalfSquared(1, 1); + // Now T and mm have been recalculated in natural units using BLHalfSquared(1,1). + + /* + double sol = + (double)(MFactor * BLHalfSquared(Bmag, Lcoh) * (1 + exp(-GammaL) - 2 * exp(-GammaL / 2) * cos(phi))); + */ + + RESTDebug << "Axion-photon transmission probability : " << sol << RESTendl; + + return (Double_t)sol; +#endif } /////////////////////////////////////////////// @@ -199,9 +328,15 @@ Double_t TRestAxionPhotonConversion::GammaTransmissionProbability(Double_t Bmag, /// /// The returned value is given for g_ag = 10^-10 GeV-1 /// -Double_t TRestAxionPhotonConversion::AxionAbsorptionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, - Double_t ma, Double_t mg, - Double_t absLength) { +Double_t TRestAxionField::AxionAbsorptionProbability(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, + Double_t mg, Double_t absLength) { +#ifndef USE_MPFR + RESTWarning + << "MPFR libraries not linked to REST libraries. Try adding -DREST_MPFr=ON to your REST compilation" + << RESTendl; + RESTWarning << "TRestAxionField::GammaTransmissionProbability will return 0" << RESTendl; + return 0; +#else mpfr::mpreal axionMass = ma; mpfr::mpreal cohLength = Lcoh / 1000.; // Default REST units are mm; @@ -211,8 +346,7 @@ Double_t TRestAxionPhotonConversion::AxionAbsorptionProbability(Double_t Bmag, D if (fDebug) { RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " TRestAxionPhotonConversion::GammaTransmissionProbability. Parameter summary" - << RESTendl; + RESTDebug << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl; RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; RESTDebug << " Axion mass : " << ma << " eV" << RESTendl; RESTDebug << " Axion energy : " << Ea << " keV" << RESTendl; @@ -225,7 +359,7 @@ Double_t TRestAxionPhotonConversion::AxionAbsorptionProbability(Double_t Bmag, D if (ma == 0.0 && photonMass == 0.0) return BLHalfSquared(Bmag, Lcoh); mpfr::mpreal q = (ma * ma - photonMass * photonMass) / 2. / Ea / 1000.0; - mpfr::mpreal l = cohLength * PhMeterIneV; + mpfr::mpreal l = cohLength * REST_Physics::PhMeterIneV; mpfr::mpreal phi = q * l; mpfr::mpreal Gamma = absLength; @@ -258,11 +392,12 @@ Double_t TRestAxionPhotonConversion::AxionAbsorptionProbability(Double_t Bmag, D if (fDebug) RESTDebug << "Axion-photon absorption probability : " << sol << RESTendl; return sol; +#endif } /// Commented because it uses ComplexReal structure that is moved to TRestAxionFieldPropagationProcess class /* -void TRestAxionPhotonConversion::PropagateAxion(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, +void TRestAxionField::PropagateAxion(Double_t Bmag, Double_t Lcoh, Double_t Ea, Double_t ma, Double_t mg, Double_t absLength) { mpfr::mpreal axionMass = ma; mpfr::mpreal cohLength = Lcoh / 1000.; // Default REST units are mm; @@ -272,7 +407,7 @@ void TRestAxionPhotonConversion::PropagateAxion(Double_t Bmag, Double_t Lcoh, Do if (fDebug) { RESTDebug << "+--------------------------------------------------------------------------+" << -RESTendl; RESTDebug << " TRestAxionPhotonConversion::GammaTransmissionProbability. Parameter summary" << +RESTendl; RESTDebug << " TRestAxionField::GammaTransmissionProbability. Parameter summary" << RESTendl; RESTDebug << " Photon mass : " << photonMass << " eV" << RESTendl; RESTDebug << " Axion mass : " << ma << " eV" << RESTendl; RESTDebug << " Axion energy : " << Ea << " keV" << RESTendl; RESTDebug << " Lcoh : " << Lcoh << " mm" << RESTendl; RESTDebug << " Bmag : " << Bmag << " T" << RESTendl; RESTDebug << diff --git a/src/TRestAxionFieldPropagationProcess.cxx b/src/TRestAxionFieldPropagationProcess.cxx index 96b6e4aa..afc6acf3 100644 --- a/src/TRestAxionFieldPropagationProcess.cxx +++ b/src/TRestAxionFieldPropagationProcess.cxx @@ -50,6 +50,43 @@ /// In a first approach this process will be only valid for the axion propagation inside a single magnetic /// volume, until it is confirmed the process is valid for any number of volumes. /// +/// This process requires at least the definition of a magnetic field using a TRestAxionMagneticField +/// metadata definition. Optionally, if the field uses a buffer gas we may define a TRestAxionBufferGas +/// that will include all the necessary gas properties, such as photon mass or absorption. +/// +/// Two metadata parameters can be defined inside this process: +/// +/// * **integrationStep** (Default:50mm): The integration length used for field integration along the particle +/// track. +/// * **bufferGasAdditionalLength** (Default:0mm): In the case we use a buffer gas we may include an +/// additional length that the particle needs to travel inside that buffer gas but outside the magnetic field +/// volume, the result will be written as an independent efficiency at the `transmission` observable inside +/// the analysis tree. +/// +/// List of observables generated by this process: +/// +/// * **fieldAverage**: The averaged magnetic field calculated along the particle track. +/// * **probability**: The final axion-photon conversion probability. +/// * **coherenceLength**: The lentgh of magnetic field region traversed by the particle. +/// * **transmission**: This observable will register the photon transmission produced by an additional +/// buffer gas length at the end of the magnetic region. +/// +/// This process can be tested using the RML files found inside the directory +/// `/pipeline/ray-tracing/axion-field/`. The following commands will generate plots with different +/// distributions containing different axion-photon probability calculations. +/// +/// \code +/// restManager --c photonConversion.rml +/// restManager --c plots.rml --f AxionPhotonProbability.root +/// \endcode +/// +/// The previous commands were used to generate the following plots where a total of 100,000 events +/// were generated. +/// +/// \htmlonly \endhtmlonly +/// +/// ![Ray-tracing distributions of probability, field average and coherence length](axionFieldPlots.png) +/// ///-------------------------------------------------------------------------- /// /// RESTsoft - Software for Rare Event Searches with TPCs @@ -65,6 +102,9 @@ /// 2020-March: Review and validation of this process. /// Javier Galan and Krešimir Jakovčić /// +/// 2022-November: Finally including non-homogeneous field integration +/// Javier Galan +/// /// /// \class TRestAxionFieldPropagationProcess /// \author Javier Galan @@ -75,6 +115,7 @@ #include "TRestAxionFieldPropagationProcess.h" #include +#include using namespace std; ClassImp(TRestAxionFieldPropagationProcess); @@ -117,38 +158,70 @@ void TRestAxionFieldPropagationProcess::Initialize() { /// should be initialized here. /// void TRestAxionFieldPropagationProcess::InitProcess() { - RESTDebug << "Entering ... TRestAxionGeneratorProcess::InitProcess" << RESTendl; + RESTDebug << "Entering ... TRestAxionFieldPropagationProcess::InitProcess" << RESTendl; - fField = (TRestAxionMagneticField*)this->GetMetadata("TRestAxionMagneticField"); + fMagneticField = (TRestAxionMagneticField*)this->GetMetadata("TRestAxionMagneticField"); - if (!fField) { + RESTDebug << "Magnetic field : " << fMagneticField << RESTendl; + + if (!fMagneticField) { RESTError << "TRestAxionFieldPropagationprocess. Magnetic Field was not defined!" << RESTendl; exit(0); } + + if (!fAxionField) { + fAxionField = new TRestAxionField(); + + fBufferGas = (TRestAxionBufferGas*)this->GetMetadata("TRestAxionBufferGas"); + if (fBufferGas) fAxionField->AssignBufferGas(fBufferGas); + } + + RESTDebug << "Axion-field : " << fAxionField << RESTendl; + RESTDebug << "Buffer gas : " << fBufferGas << RESTendl; } TRestEvent* TRestAxionFieldPropagationProcess::ProcessEvent(TRestEvent* evInput) { - // Already done by TRestAxionEventProcess fAxionEvent = (TRestAxionEvent*)evInput; + RESTDebug << "TRestAxionFieldPropagationProcess::ProcessEvent : " << fAxionEvent->GetID() << RESTendl; std::vector trackBounds = - fField->GetFieldBoundaries(fAxionEvent->GetPosition(), fAxionEvent->GetDirection()); + fMagneticField->GetFieldBoundaries(fAxionEvent->GetPosition(), fAxionEvent->GetDirection()); - if (trackBounds.size() == 0) { - fAxionEvent->SetBSquared(0); - fAxionEvent->SetLConversion(0); - } else { - Double_t B = fField->GetTransversalFieldAverage(trackBounds[0], trackBounds[1], 100, 200); - Double_t lenght = (trackBounds[1] - trackBounds[0]).Mag(); + Double_t prob = 0; + Double_t lCoh = 0; + Double_t transmission = 0; + Double_t fieldAverage = 0; + if (trackBounds.size() == 2) { + std::vector bProfile = fMagneticField->GetTransversalComponentAlongPath( + trackBounds[0], trackBounds[1], fIntegrationStep); - fAxionEvent->SetBSquared(B * B); - fAxionEvent->SetLConversion(lenght); + fieldAverage = std::accumulate(bProfile.begin(), bProfile.end(), 0.0) / bProfile.size(); + + Double_t Ea = fAxionEvent->GetEnergy(); + Double_t ma = fAxionEvent->GetMass(); + + prob = fAxionField->GammaTransmissionProbability(bProfile, fIntegrationStep, Ea, ma); + + lCoh = (bProfile.size() - 1) * fIntegrationStep; + + if (fBufferGas && fBufferGasAdditionalLength > 0) { + Double_t Gamma = fBufferGas->GetPhotonAbsorptionLength(Ea); // cm-1 + Double_t GammaL = Gamma * lCoh * units("cm"); + transmission = exp(-GammaL); + } } + SetObservableValue("fieldAverage", fieldAverage); + SetObservableValue("probability", prob); + SetObservableValue("coherenceLength", lCoh); + SetObservableValue("transmission", transmission); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fAxionEvent->PrintEvent(); /// Missing to propagate the axion to the end of magnet bore? + /// May not be necessary, it can be done by TRestAxionTransportProcess if user needs any process + /// that changes direction, this is done for example by optics processes internally return fAxionEvent; } diff --git a/src/TRestAxionLikelihood.cxx b/src/TRestAxionLikelihood.cxx index dac3d9c4..f5d26c17 100644 --- a/src/TRestAxionLikelihood.cxx +++ b/src/TRestAxionLikelihood.cxx @@ -84,9 +84,9 @@ void TRestAxionLikelihood::Initialize() { fBufferGas = new TRestAxionBufferGas(); // Conversion probabilities as defined by van Bibber paper - fPhotonConversion = new TRestAxionPhotonConversion(); + fAxionField = new TRestAxionField(); - fPhotonConversion->AssignBufferGas(fBufferGas); + fAxionField->AssignBufferGas(fBufferGas); // Solar axion flux on earth fAxionSpectrum = new TRestAxionSpectrum(fConfigFileName.c_str()); @@ -288,8 +288,8 @@ Double_t TRestAxionLikelihood::GetSignal(Double_t ma, Double_t g10_4, Double_t r Double_t dE = 0.01; for (Double_t en = fErange.X(); en < fErange.Y(); en = en + dE) { Double_t Phi_a = fAxionSpectrum->GetDifferentialSolarAxionFlux(en); - // TODO Needs to be readapted to the new TRestAxionPhotonConversion implementation - // Double_t Pa_g = fPhotonConversion->GammaTransmissionProbability(en, fBmag, ma); + // TODO Needs to be readapted to the new TRestAxionField implementation + // Double_t Pa_g = fField->GammaTransmissionProbability(en, fBmag, ma); Double_t Pa_g = 0; // Just dummy 0 probability!! Double_t nGamma = Pa_g * Phi_a; @@ -347,7 +347,7 @@ Double_t TRestAxionLikelihood::LogLikelihood(Double_t ma, Double_t g10_4, Double // Double_t Ngamma = solarFluxG10 * magnetArea * (fTExpVacuum * 3600.) * fEfficiency; // Double_t nGammaVacuum = - // fPhotonConversion->SetAxionMass( Double_t m ) { fAxionMass = m; } + // fField->SetAxionMass( Double_t m ) { fAxionMass = m; } return lhood; } @@ -361,9 +361,9 @@ void TRestAxionLikelihood::InitFromConfigFile() { fRmag = GetDblParameterWithUnits("Rmag", 300); fLmag = GetDblParameterWithUnits("Lmag", 4500.); - /// This needs to be reviewed. TRestAxionPhotonConversion does not store those parameter anymore. - // fPhotonConversion->SetCoherenceLength(fLmag); - // fPhotonConversion->SetMagneticField(fBmag); + /// This needs to be reviewed. TRestAxionField does not store those parameter anymore. + // fField->SetCoherenceLength(fLmag); + // fField->SetMagneticField(fBmag); fEfficiency = StringToDouble(GetParameter("efficiency", "1")); diff --git a/src/TRestAxionOptics.cxx b/src/TRestAxionOptics.cxx index 40382aed..187560f7 100644 --- a/src/TRestAxionOptics.cxx +++ b/src/TRestAxionOptics.cxx @@ -261,7 +261,9 @@ Int_t TRestAxionOptics::TransportToEntrance(const TVector3& pos, const TVector3& fEntrancePosition = REST_Physics::MoveToPlane(pos, dir, TVector3(0, 0, 1), TVector3(0, 0, GetEntrancePositionZ())); - return fEntranceMask->GetRegion(fEntrancePosition.X(), fEntrancePosition.Y()); + Double_t x = fEntrancePosition.X(); + Double_t y = fEntrancePosition.Y(); + return fEntranceMask->GetRegion(x, y); } /////////////////////////////////////////////// @@ -290,7 +292,9 @@ Int_t TRestAxionOptics::TransportToMiddle(const TVector3& pos, const TVector3& d fMiddlePosition = REST_Physics::MoveToPlane(pos, dir, TVector3(0, 0, 1), TVector3(0, 0, 0)); fMiddleDirection = dir; - return fMiddleMask->GetRegion(fMiddlePosition.X(), fMiddlePosition.Y()); + Double_t x = fMiddlePosition.X(); + Double_t y = fMiddlePosition.Y(); + return fMiddleMask->GetRegion(x, y); } /////////////////////////////////////////////// @@ -319,7 +323,9 @@ Int_t TRestAxionOptics::TransportToExit(const TVector3& pos, const TVector3& dir REST_Physics::MoveToPlane(pos, dir, TVector3(0, 0, 1), TVector3(0, 0, GetExitPositionZ())); fExitDirection = dir; - return fExitMask->GetRegion(fExitPosition.X(), fExitPosition.Y()); + Double_t x = fExitPosition.X(); + Double_t y = fExitPosition.Y(); + return fExitMask->GetRegion(x, y); } /////////////////////////////////////////////// diff --git a/src/TRestAxionOpticsMirror.cxx b/src/TRestAxionOpticsMirror.cxx index b86d0b86..e75f2ce3 100644 --- a/src/TRestAxionOpticsMirror.cxx +++ b/src/TRestAxionOpticsMirror.cxx @@ -29,7 +29,8 @@ /// The following metadata parameters are used to define the mirror /// properties: /// * **mirrorType**: It defines the mirror type (Single, Thick, Bilayer, -/// Multilayer). For the moment only `Single` layer type has been implemented. +/// Multilayer). For the moment only `Single` and `Bilayer` types have +/// been implemented. /// * **layer**: It defines the layer material. Chemical formula. /// * **layerThickness**: It defines the layer thickness in nm. /// * **substrate**: It defines the substrate material. Chemical formula. diff --git a/src/TRestAxionOpticsProcess.cxx b/src/TRestAxionOpticsProcess.cxx index 8e418d33..85df0c0d 100644 --- a/src/TRestAxionOpticsProcess.cxx +++ b/src/TRestAxionOpticsProcess.cxx @@ -115,6 +115,17 @@ void TRestAxionOpticsProcess::InitProcess() { } else { RESTError << "TRestAxionOptics::InitProcess. No sucess instantiating optics." << RESTendl; } + + if (fOpticalAxis) { + // The outgoing particle will be referenced to the optical axis. + // The particle will not be counter-rotated when the process ends + SkipEndProcessRotation(); + } else { + // The outgoing particle will be referenced to the universal axis + // The optics will apply a deflection angle, but the main axis will be kept aligned with the previous + // axis reference + SkipEndProcessRotation(false); + } } /////////////////////////////////////////////// @@ -148,16 +159,3 @@ TRestEvent* TRestAxionOpticsProcess::ProcessEvent(TRestEvent* evInput) { return fAxionEvent; } - -////////////////////////////////////////////////////////////////////////// -/// \brief End of event process. -/// -void TRestAxionOpticsProcess::EndOfEventProcess(TRestEvent* evInput) { - if (fOpticalAxis) { - // The outgoing particle will be referenced to the optical axis - TRestAxionEventProcess::EndOfEventProcess(evInput); - } else { - // The outgoing particle will be referenced to the universal axis - TRestEventProcess::EndOfEventProcess(evInput); - } -} diff --git a/src/TRestAxionTransportProcess.cxx b/src/TRestAxionTransportProcess.cxx index ef901b19..4894228b 100644 --- a/src/TRestAxionTransportProcess.cxx +++ b/src/TRestAxionTransportProcess.cxx @@ -125,9 +125,8 @@ TRestEvent* TRestAxionTransportProcess::ProcessEvent(TRestEvent* evInput) { TVector3 inDir = fAxionEvent->GetDirection(); if ((inPos.Z() - fZPosition) * inDir.Z() >= 0) { - RESTWarning - << "TRestAxionTransportProcess::ProcessEvent. Not appropiate particle direction to reach Z=" - << fZPosition << RESTendl; + RESTWarning << "TRestAxionTransportProcess::ProcessEvent (" << GetName() + << "). Not appropiate particle direction to reach Z=" << fZPosition << RESTendl; RESTWarning << "Axion position. Z:" << inPos.Z() << " direction Z: " << inDir.Z() << RESTendl; fAxionEvent->PrintEvent();