diff --git a/.github/workflows/validation.yml b/.github/workflows/validation.yml index ef47fa82..ee2180b2 100644 --- a/.github/workflows/validation.yml +++ b/.github/workflows/validation.yml @@ -80,7 +80,7 @@ jobs: restRoot -b -q GetMagneticField_test.C cd ../boundary/ restRoot -b -q Boundaries_test.C - - name: Optics + - name: Optics Mirrors run: | source ${{ env.REST_PATH }}/thisREST.sh cd pipeline/metadata/optics/ @@ -91,8 +91,12 @@ jobs: wget https://sultan.unizar.es/axionlib-data/opticsMirror/Reflectivity_Single_Au_250_Ni_0.4.N901f wget https://sultan.unizar.es/axionlib-data/opticsMirror/Transmission_Single_Au_250_Ni_0.4.N901f python mirrors.py - #python optics.py - #python basic.py + - name: MCPL Optics + run: | + source ${{ env.REST_PATH }}/thisREST.sh + cd pipeline/metadata/optics/ + wget https://sultan.unizar.es/axionlib-data/optics/optics.rml + python mcpl.py - name: X-ray transmission run: | source ${{ env.REST_PATH }}/thisREST.sh @@ -114,3 +118,29 @@ jobs: python solarTests.py python solarPlot.py python compare.py + + Ray-tracing: + name: Check ray-tracing processing chain + 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: XMM optics bench ray-tracing + run: | + source ${{ env.REST_PATH }}/thisREST.sh + export REST_NEVENTS=1000 + export REST_RUN=100 + wget https://sultan.unizar.es/axionlib-data/optics/xmm.rml + wget https://sultan.unizar.es/axionlib-data/optics/XMM.Wolter + wget https://sultan.unizar.es/axionlib-data/opticsMirror/Reflectivity_Single_Au_250_Ni_0.4.N901f + 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")' diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 56b1f97b..fce91ddb 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -5,6 +5,7 @@ stages: - build - loadRESTLibs - metadata + - ray-tracing before_script: - export USER="axion" @@ -104,6 +105,17 @@ mirrors: variables: - $CRONJOB +MCPL Optics: + stage: metadata + script: + - . ${CI_PROJECT_DIR}/install/thisREST.sh + - cd ${CI_PROJECT_DIR}/pipeline/metadata/optics/ + - wget https://sultan.unizar.es/axionlib-data/optics/optics.rml + - python mcpl.py + except: + variables: + - $CRONJOB + xRayTransmission: stage: metadata script: @@ -131,3 +143,20 @@ solarFlux: except: variables: - $CRONJOB + +XMM optics bench ray-tracing: + stage: ray-tracing + script: + - . ${CI_PROJECT_DIR}/install/thisREST.sh + - cd ${CI_PROJECT_DIR}/pipeline/ray-tracing/ + - export REST_NEVENTS=1000 + - export REST_RUN=100 + - wget https://sultan.unizar.es/axionlib-data/opticsMirror/Reflectivity_Single_Au_250_Ni_0.4.N901f + - wget https://sultan.unizar.es/axionlib-data/opticsMirror/Transmission_Single_Au_250_Ni_0.4.N901f + - wget https://sultan.unizar.es/axionlib-data/optics/xmm.rml + - wget https://sultan.unizar.es/axionlib-data/optics/XMM.Wolter + - restManager --c opticsBench.rml + - restRoot ValidateXMM.C'("OpticsBench_Yaw_0.05_Dev_0.005_BabyIAXO_Run00100.root")' + except: + variables: + - $CRONJOB diff --git a/CMakeLists.txt b/CMakeLists.txt index 120d1617..27dd9055 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -68,5 +68,10 @@ INSTALL(DIRECTORY ./data/ COMPONENT install ) +INSTALL(DIRECTORY ./examples/ + DESTINATION ./examples/axion/ + COMPONENT install + ) + file(GLOB_RECURSE MAC "${CMAKE_CURRENT_SOURCE_DIR}/macros/*") INSTALL(FILES ${MAC} DESTINATION ./macros/axion) diff --git a/README.md b/README.md index be1655a3..3cc4f3ee 100644 --- a/README.md +++ b/README.md @@ -59,3 +59,10 @@ cd framework/build cmake -DRESTLIB_AXION=ON ../ make -j4 install ``` + +### Publications + +This repository makes use of the following published codes: +- K. Altenmuller et al, REST-for-Physics, a ROOT-based framework for event oriented data analysis and combined Monte Carlo response, [Computer Physics Communications 273, April 2022, 108281](https://doi.org/10.1016/j.cpc.2021.108281). +- S.Hoof, J.Jaeckel, T.J.Lennert, Quantifying uncertainties in the solar axion flux and their impact on determining axion model parameters, [JCAP09(2021)006](https://doi.org/10.1088/1475-7516/2021/09/006). +- T.Kittelmann, E.Klinkby, E.B.Knudsen, P.Willendrup, X.X.Cai, K.Kanaki, Monte Carlo Particle Lists: MCPL, Computer Physics Communications 218 (2017) 17–42](http://dx.doi.org/10.17632/cby92vsv5g.1). diff --git a/data b/data index f3bc577c..c7ad6a5c 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit f3bc577c76f2c5af6fcc1a30f9f0288b6887cfd9 +Subproject commit c7ad6a5cf338989786532b76403c63183321550a diff --git a/examples/README.md b/examples/README.md index 78424d33..d39aef86 100644 --- a/examples/README.md +++ b/examples/README.md @@ -1 +1,3 @@ TOBE written, say few words listing, and describing the examples, and briefly explain how to use them. + +- **opticsBench.rml**: It defines a simple generator with a parallel flux in order to validate the optics process. diff --git a/examples/axionGenerator.rml b/examples/axionGenerator.rml index 037c81cf..d8088f69 100644 --- a/examples/axionGenerator.rml +++ b/examples/axionGenerator.rml @@ -1,113 +1,80 @@ - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - diff --git a/examples/opticsBench.rml b/examples/opticsBench.rml new file mode 100644 index 00000000..f839b1a0 --- /dev/null +++ b/examples/opticsBench.rml @@ -0,0 +1,70 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/opticsPlots.rml b/examples/opticsPlots.rml new file mode 100644 index 00000000..ca55e8b7 --- /dev/null +++ b/examples/opticsPlots.rml @@ -0,0 +1,80 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inc/TRestAxionAnalysisProcess.h b/inc/TRestAxionAnalysisProcess.h index 13c18cea..41c70809 100644 --- a/inc/TRestAxionAnalysisProcess.h +++ b/inc/TRestAxionAnalysisProcess.h @@ -32,8 +32,6 @@ class TRestAxionAnalysisProcess : public TRestEventProcess { /// A pointer to the specific TRestAxionEvent TRestAxionEvent* fAxionEvent; //! - void InitFromConfigFile() override; - void Initialize() override; void LoadDefaultConfig(); diff --git a/inc/TRestAxionDeviationProcess.h b/inc/TRestAxionDeviationProcess.h new file mode 100644 index 00000000..b66355bb --- /dev/null +++ b/inc/TRestAxionDeviationProcess.h @@ -0,0 +1,76 @@ +/************************************************************************* + * 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 RestCore_TRestAxionDeviationProcess +#define RestCore_TRestAxionDeviationProcess + +#include "TRandom3.h" + +#include "TRestAxionEvent.h" +#include "TRestAxionEventProcess.h" + +//! A process to deviate the axion direction by a given yaw and pitch angle distributions +class TRestAxionDeviationProcess : public TRestAxionEventProcess { + private: + /// The angle that defines the cone directrix that defines the maximum deviation + Double_t fDevAngle = 0; //< + + /// Seed used in random generator + Int_t fSeed = 0; //< + + /// Internal process random generator + TRandom3* fRandom = nullptr; //! + + void LoadDefaultConfig(); + + protected: + public: + void InitProcess() override; + + void Initialize() override; + + TRestEvent* ProcessEvent(TRestEvent* evInput) override; + + void LoadConfig(std::string cfgFilename, std::string name = ""); + + /// It prints out the process parameters stored in the metadata structure + void PrintMetadata() override { + BeginPrintProcess(); + + RESTMetadata << "Cone directrix angle: " << fDevAngle * units("degrees") << " degrees" << RESTendl; + + EndPrintProcess(); + } + + /// Returns the name of this process + const char* GetProcessName() const override { return "axionDeviation"; } + + // Constructor + TRestAxionDeviationProcess(); + TRestAxionDeviationProcess(char* cfgFileName); + + // Destructor + ~TRestAxionDeviationProcess(); + + ClassDefOverride(TRestAxionDeviationProcess, 1); +}; +#endif diff --git a/inc/TRestAxionEvent.h b/inc/TRestAxionEvent.h index 2bbc1194..183f95da 100644 --- a/inc/TRestAxionEvent.h +++ b/inc/TRestAxionEvent.h @@ -30,51 +30,78 @@ #include "TVector3.h" #include "TRestEvent.h" +#include "TRestSystemOfUnits.h" /// An event data class to define the parameters related to an axion particle class TRestAxionEvent : public TRestEvent { private: - TVector3 fPosition; //-> Particle position - TVector3 fDirection; //-> Unitary direction of movement - Double_t fEnergy = 0; //-> Energy of axion in keV + /// Particle position + TVector3 fPosition; //< - Double_t fMass = 0.; //-> Axion mass in eV + /// Particle momentum. Unitary direction. + TVector3 fDirection; //< - Double_t fGammaProbability = 0; //-> The conversion probability P_{ag} + /// Initial energy of the axion. + Double_t fEnergy = 0; //< - Double_t fEfficiency = 1; //-> To include any loss of signal transmission/efficiency + /// 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; //< protected: public: - TVector3* GetPosition() { return &fPosition; } + TVector3 GetPosition() { return fPosition; } Double_t GetPositionX() { return fPosition.X(); } // returns value in mm Double_t GetPositionY() { return fPosition.Y(); } // returns value in mm Double_t GetPositionZ() { return fPosition.Z(); } // returns value in mm - TVector3* GetDirection() { return &fDirection; } + TVector3 GetDirection() { return fDirection; } Double_t GetDirectionX() { return fDirection.X(); } // returns normalized vector x-component. Double_t GetDirectionY() { return fDirection.Y(); } // returns normalized vector y-component Double_t GetDirectionZ() { return fDirection.Z(); } // returns normalized vector z-component - Double_t GetEnergy() { return fEnergy; } // returns value in keV - Double_t GetMass() { return fMass; } // returns value in eV - Double_t GetEfficiency() { return fEfficiency; } + 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 GetGammaProbability() { return fGammaProbability; } + 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)); } + void Translate(const TVector3& delta); + void SetDirection(TVector3 dir) { fDirection = dir; } void SetDirection(Double_t px, Double_t py, Double_t pz) { SetDirection(TVector3(px, py, pz)); } + void RotateZX(const TVector3& center, Double_t phi, Double_t theta); + void RotateXZ(const TVector3& center, Double_t theta, Double_t phi); + + void RotateXY(const TVector3& center, Double_t pitch, Double_t yaw); + void RotateYX(const TVector3& center, Double_t yaw, Double_t pitch); + void SetEnergy(Double_t en) { fEnergy = en; } void SetMass(Double_t m) { fMass = m; } - void SetGammaProbability(Double_t p) { fGammaProbability = p; } - void SetEfficiency(Double_t eff) { fEfficiency = eff; } + 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; } virtual void Initialize(); @@ -82,9 +109,7 @@ class TRestAxionEvent : public TRestEvent { TPad* DrawEvent(TString option = ""); - // Construtor TRestAxionEvent(); - // Destructor ~TRestAxionEvent(); ClassDef(TRestAxionEvent, 1); diff --git a/inc/TRestAxionEventProcess.h b/inc/TRestAxionEventProcess.h new file mode 100644 index 00000000..df6147c3 --- /dev/null +++ b/inc/TRestAxionEventProcess.h @@ -0,0 +1,67 @@ +/************************************************************************* + * 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 RestCore_TRestAxionEventProcess +#define RestCore_TRestAxionEventProcess + +#include "TRestAxionEvent.h" +#include "TRestEventProcess.h" + +/// A base class for any axion event process. Defines position, rotation and component displacement. +class TRestAxionEventProcess : public TRestEventProcess { + private: + /// The position respect which the rotation will be applied + TVector3 fCenter = TVector3(0, 0, 0); + + /// The rotation angle respect to the Y-axis + Double_t fYaw = 0; + + /// The rotation angle with respect to X-axis + Double_t fPitch = 0; + + protected: + /// A pointer to the specific TRestAxionEvent + TRestAxionEvent* fAxionEvent; //! + + void BeginPrintProcess(); + void EndPrintProcess(); + + TVector3 GetCenter() const { return fCenter; } + + public: + RESTValue GetInputEvent() const override { return fAxionEvent; } + RESTValue GetOutputEvent() const override { return fAxionEvent; } + + virtual void InitProcess() override {} + + /// Begin of event process, preparation work. Called right before ProcessEvent() + virtual void BeginOfEventProcess(TRestEvent* evInput = nullptr) override; + + /// End of event process. Called directly after ProcessEvent() + virtual void EndOfEventProcess(TRestEvent* evInput = nullptr) override; + + TRestAxionEventProcess(); + ~TRestAxionEventProcess(); + + ClassDefOverride(TRestAxionEventProcess, 1); +}; +#endif diff --git a/inc/TRestAxionFieldPropagationProcess.h b/inc/TRestAxionFieldPropagationProcess.h index d61e796f..b617705a 100644 --- a/inc/TRestAxionFieldPropagationProcess.h +++ b/inc/TRestAxionFieldPropagationProcess.h @@ -27,216 +27,25 @@ #include "TVectorD.h" #include "TRestAxionEvent.h" +#include "TRestAxionEventProcess.h" #include "TRestAxionMagneticField.h" -#include "TRestAxionPhotonConversion.h" -#include "TRestEventProcess.h" #include "TRestPhysics.h" -#include "mpreal.h" -/// A structure to define the two components of a complex number using real precision. -/// To be used inside TRestAxionFieldPropagationProcess. Copied from TRestAxionPhotonConversion.h -struct ComplexReal { - /// The real part of the number - mpfr::mpreal real = 0; - - /// The imaginary part of the number - mpfr::mpreal img = 0; -}; - -//! A process to introduce the axion-photon conversion probability in the signal generation chain -class TRestAxionFieldPropagationProcess : public TRestEventProcess { +//! A process to introduce the magnetic field profile integration along the track +class TRestAxionFieldPropagationProcess : public TRestAxionEventProcess { private: - /// complex axion field amplitude. - ComplexReal faxionAmplitude; //! - - /// complex amplitude of the photon field component parallel to the transverse magnetic field. - ComplexReal fparallelPhotonAmplitude; //! - - /// complex amplitude of the photon field component orthogonal to the transverse magnetic field. - ComplexReal forthogonalPhotonAmplitude; //! - - ///////////////////////////////////////////////////////////////////////// - // ----- 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; - } - - void PrintComplex(ComplexReal); - - void InitFromConfigFile() override; + /// A pointer to the magnetic field description stored in TRestRun + TRestAxionMagneticField* fField; //! void Initialize() override; - void LoadDefaultConfig(); - - /// A pointer to the specific TRestAxionEvent - TRestAxionEvent* fAxionEvent; //! - - /// A pointer to the magnetic field stored in TRestRun - TRestAxionMagneticField* fAxionMagneticField; //! - - /// A pointer for the gamma conversion probability - TRestAxionPhotonConversion* fAxionPhotonConversion; //! - - /// A pointer for the gamma conversion probability - TRestAxionBufferGas* fAxionBufferGas; //! - - /// Variables for the output AxionEvent position - TString fMode; - TVector3 fFinalNormalPlan; - TVector3 fFinalPositionPlan; - Double_t fDistance; - - protected: public: void InitProcess() override; - any GetInputEvent() const override { return fAxionEvent; } - any GetOutputEvent() const override { return fAxionEvent; } - - TRestEvent* ProcessEvent(TRestEvent* evInput) override; - - void LoadConfig(std::string cfgFilename, std::string name = ""); - - /// It prints out the process parameters stored in the metadata structure - void PrintMetadata() override { - BeginPrintProcess(); - - RESTMetadata << "mode: " << fMode << RESTendl; - - Double_t x = fFinalPositionPlan.X(); - Double_t y = fFinalPositionPlan.Y(); - Double_t z = fFinalPositionPlan.Z(); - RESTMetadata << "finalPositionPlan = ( " << x << " ," << y << ", " << z << ") mm" << RESTendl; - - x = fFinalNormalPlan.X(); - y = fFinalNormalPlan.Y(); - z = fFinalNormalPlan.Z(); - RESTMetadata << "finalNormalPlan = ( " << x << ", " << y << ", " << z << ")" << RESTendl; - - RESTMetadata << "Distance: " << fDistance << " mm" << RESTendl; - - EndPrintProcess(); - } - - bool IsInBoundedPlan(TVector3 pos, Int_t i, Int_t p); - std::vector InOut(std::vector bounds, TVector3 dir); - std::vector FindBoundariesOneVolume(TVector3 pos, TVector3 dir, Int_t p); - std::vector FieldBoundary(std::vector boundaries, Double_t minStep); - - /// Returns the boundaries of the axion passed through magnetic fields - std::vector> FindFieldBoundaries(Double_t minStep = -1); - - // TVectorD GetFieldVector(TVector3 in, TVector3 out, Int_t N = 0); - - // Obsolete - // TVector3 MoveToFinalDistance(TVector3 pos, TVector3 dir, Double_t distance); - // TVector3 MoveToPlan(TVector3 pos, TVector3 dir, Double_t f, Int_t i); - // TVector3 FinalPositionInPlan(TVector3 pos, TVector3 dir, TVector3 normalPlan, TVector3 pointPlan); - - // Calculates amplitudes of the axion field, parallel component of the photon field and orthogonal - // component of the photon field (with respect to the transverse component of the magnetic field) - // after passing one segment of the particle trajectory along which it is assumed - // that the transverse component of the magnetic field with respect to the particle propagation direction - // is not equal zero. - void CalculateAmplitudesInSegment(ComplexReal& faxionAmplitude, ComplexReal& fparallelPhotonAmplitude, - ComplexReal& forthogonalPhotonAmplitude, TVector3& averageBT, - mpfr::mpreal axionMass, mpfr::mpreal photonMass, mpfr::mpreal Ea, - TVector3 from, TVector3 to, mpfr::mpreal CommonPhase, - mpfr::mpreal OrthogonalPhase); - - // Calculates amplitudes of the axion field, parallel component of the photon field and orthogonal - // component - // of the photon field after passing one subsegment of the particle trajectory along which it is assumed - // that the magnetic field is constant - void CalculateAmplitudesInSubsegment(ComplexReal& axionAmplitude, ComplexReal& parallelPhotonAmplitude, - ComplexReal& orhogonalPhotonAmplitude, mpfr::mpreal theta, - mpfr::mpreal lambda, Double_t length, mpfr::mpreal CommonPhase, - mpfr::mpreal OrthogonalPhase); + RESTValue GetInputEvent() const override { return fAxionEvent; } + RESTValue GetOutputEvent() const override { return fAxionEvent; } - // Calculates amplitudes of the axion field, parallel component of the photon field and orthogonal - // component - // of the photon field after passing one segment of the particle trajectory along which the transverse - // component of the magnetic field is zero - void PropagateWithoutBField(ComplexReal& axionAmplitude, ComplexReal& parallelPhotonAmplitude, - ComplexReal& orthogonalPhotonAmplitude, mpfr::mpreal axionMass, - mpfr::mpreal photonMass, mpfr::mpreal Ea, TVector3 from, TVector3 to); + TRestEvent* ProcessEvent(TRestEvent* eventInput) override; /// Returns the name of this process const char* GetProcessName() const override { return "axionFieldPropagation"; } diff --git a/inc/TRestAxionGeneratorProcess.h b/inc/TRestAxionGeneratorProcess.h index 44974746..392cde52 100644 --- a/inc/TRestAxionGeneratorProcess.h +++ b/inc/TRestAxionGeneratorProcess.h @@ -26,70 +26,47 @@ #include "TRestAxionEvent.h" #include "TRestEventProcess.h" -#include "TRestAxionSpectrum.h" +#include "TRestAxionSolarFlux.h" #include "TRandom3.h" -//! A process to generate axions following a particular solar axion model +//! A process to initialize the axion event (mainly through TRestAxionSolarFlux) class TRestAxionGeneratorProcess : public TRestEventProcess { private: /// A pointer to the specific TRestAxionEvent output TRestAxionEvent* fOutputAxionEvent; //! + /// A pointer to the TRestAxionSolarFlux metadata description + TRestAxionSolarFlux* fAxionFlux; //! + /// Used internally to define the event id Int_t fCounter = 0; //! - /// A pointer to the axion model stored in TRestRun - TRestAxionSpectrum* fAxionSpectrum; //! - - /// Random number generator - TRandom3* fRandom; //! - - /// The solar energy range - TVector2 fEnergyRange; //-> - - /// Energy step - Double_t fEnergyStep; //-> - - /// Axion mass - Double_t fAxionMass; //-> + /// Internal process random generator + TRandom3* fRandom = nullptr; //! - /// The angular distribution generator type - TString fAngularDistribution; //-> + /// The axion mass + Double_t fAxionMass = 0; //< - /// The main direction of the angular distribution - TVector3 fAngularDirection; //-> + /// The target size in mm (or generator source extension) for the generator. + Double_t fTargetRadius = 800; //< - /// The spatial distribution generator type - TString fSpatialDistribution; //-> + /// The generator type (solarFlux/flat) + TString fGeneratorType = "solarFlux"; //< - /// The spatial distribution generator type - Double_t fSpatialRadius; //-> + /// Seed used in random generator + Int_t fSeed = 0; //< - /// The spatial origin from the spatial distribution - TVector3 fSpatialOrigin; //-> + /// It defines the minimum energy as a cut-off to the generator. Default is 50eV. + Double_t fMinEnergy = 0.05; //< - /// Rotation vector (according to x,y,z) for the circle wall defined from the plan (X,Y) - insead of - /// rotating magnet Volumes - TVector3 fRotation; //-> - - /// The normal vector of the wall (other way to do a rotation) - TVector3 fNormalPlan; //-> - - /// Mode of rotated circle Wall construction - TString fMode; //-> - - void InitFromConfigFile() override; + /// It defines the maximum energy as a cut-off to the generator. No limits if equals 0. + Double_t fMaxEnergy = 0; //< void Initialize() override; void LoadDefaultConfig(); - Double_t GenerateEnergy(); - TVector3 GeneratePosition(); - TVector3 GenerateDirection(); - - protected: public: void InitProcess() override; @@ -100,32 +77,7 @@ class TRestAxionGeneratorProcess : public TRestEventProcess { void LoadConfig(std::string cfgFilename, std::string name = ""); - /// It prints out the process parameters stored in the metadata structure - void PrintMetadata() override { - BeginPrintProcess(); - - RESTMetadata << "Energy distribution" << RESTendl; - RESTMetadata << "---------------------" << RESTendl; - RESTMetadata << "Energy range : (" << fEnergyRange.X() << ", " << fEnergyRange.Y() << ") keV" - << RESTendl; - RESTMetadata << "Energy step : " << fEnergyStep << " keV" << RESTendl; - RESTMetadata << " " << RESTendl; - - RESTMetadata << "Angular distribution" << RESTendl; - RESTMetadata << "----------------------" << RESTendl; - RESTMetadata << "Type : " << fAngularDistribution << RESTendl; - RESTMetadata << "Main direction : (" << fAngularDirection.X() << "," << fAngularDirection.Y() << "," - << fAngularDirection.Z() << ")" << RESTendl; - RESTMetadata << " " << RESTendl; - RESTMetadata << "Spatial distribution" << RESTendl; - RESTMetadata << "----------------------" << RESTendl; - RESTMetadata << "Type : " << fSpatialDistribution << RESTendl; - RESTMetadata << "Radius : " << fSpatialRadius << " mm" << RESTendl; - RESTMetadata << "Origin : (" << fSpatialOrigin.X() << "," << fSpatialOrigin.Y() << "," - << fSpatialOrigin.Z() << ")" << RESTendl; - - EndPrintProcess(); - } + void PrintMetadata() override; /// Returns the name of this process const char* GetProcessName() const override { return "axionGenerator"; } diff --git a/inc/TRestAxionOptics.h b/inc/TRestAxionOptics.h index 0cc08086..6289e7d8 100644 --- a/inc/TRestAxionOptics.h +++ b/inc/TRestAxionOptics.h @@ -179,6 +179,9 @@ class TRestAxionOptics : public TRestMetadata { /// Returns the exit position from the latest propagated photon TVector3 GetExitDirection() { return fExitDirection; } + TVector3 GetLastGoodPosition(); + TVector3 GetLastGoodDirection(); + /// It returns true if the photon got reflected in the first mirror Bool_t IsFirstMirrorReflection() { return fFirstInteraction; } diff --git a/inc/TRestAxionOpticsProcess.h b/inc/TRestAxionOpticsProcess.h new file mode 100644 index 00000000..e79514d3 --- /dev/null +++ b/inc/TRestAxionOpticsProcess.h @@ -0,0 +1,89 @@ +/************************************************************************* + * 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 RestCore_TRestAxionOpticsProcess +#define RestCore_TRestAxionOpticsProcess + +#include "TRestAxionEvent.h" +#include "TRestAxionEventProcess.h" +#include "TRestAxionOptics.h" + +//! A process to introduce the response from optics in the axion signal generation chain +class TRestAxionOpticsProcess : public TRestAxionEventProcess { + private: + /// A variable to determine if the new axis will be optical or universal axis + Bool_t fOpticalAxis = false; //< + + /// A pointer to the optics description defined inside TRestRun + TRestAxionOptics* fOptics; //! + + void Initialize() override; + + void LoadDefaultConfig(); + + protected: + public: + void InitProcess() override; + + /// This InitFromConfigFile could be removed as soon as PR rest-for-physics/framework#275 + /// is solved + void InitFromConfigFile() override { + TRestEventProcess::InitFromConfigFile(); + + if (ToUpper(GetParameter("opticalAxis", "false")) == "TRUE") + fOpticalAxis = true; + else + fOpticalAxis = false; + } + + 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 + void PrintMetadata() override { + BeginPrintProcess(); + + if (fOpticalAxis) + RESTMetadata << "Output particle is described respect to the optical axis" << RESTendl; + else + RESTMetadata << "Output particle is described respect to the universal axis" << RESTendl; + + EndPrintProcess(); + } + + /// Returns the name of this process + const char* GetProcessName() const override { return "axionOptics"; } + + // Constructor + TRestAxionOpticsProcess(); + TRestAxionOpticsProcess(char* cfgFileName); + + // Destructor + ~TRestAxionOpticsProcess(); + + ClassDefOverride(TRestAxionOpticsProcess, 1); +}; +#endif diff --git a/inc/TRestAxionOpticsResponseProcess.h b/inc/TRestAxionTransportProcess.h similarity index 72% rename from inc/TRestAxionOpticsResponseProcess.h rename to inc/TRestAxionTransportProcess.h index a63b3d79..1fc8af86 100644 --- a/inc/TRestAxionOpticsResponseProcess.h +++ b/inc/TRestAxionTransportProcess.h @@ -20,30 +20,27 @@ * For the list of contributors see $REST_PATH/CREDITS. * *************************************************************************/ -#ifndef RestCore_TRestAxionOpticsResponseProcess -#define RestCore_TRestAxionOpticsResponseProcess +#ifndef RestCore_TRestAxionTransportProcess +#define RestCore_TRestAxionTransportProcess #include "TRestAxionEvent.h" -#include "TRestEventProcess.h" +#include "TRestAxionEventProcess.h" -//! A process to introduce the response from optics in the axion signal generation chain -class TRestAxionOpticsResponseProcess : public TRestEventProcess { +//! A process to transport the axion to a given z-position without changing direction +class TRestAxionTransportProcess : public TRestAxionEventProcess { private: - /// A pointer to the specific TRestAxionEvent - TRestAxionEvent* fAxionEvent; //! - - void InitFromConfigFile() override; - - void Initialize() override; + /// The Z-position where we will move the particle + Double_t fZPosition; //< void LoadDefaultConfig(); protected: public: - TRestEvent* ProcessEvent(TRestEvent* evInput) override; + void InitProcess() override; - any GetInputEvent() { return fAxionEvent; } - any GetOutputEvent() { return fAxionEvent; } + void Initialize() override; + + TRestEvent* ProcessEvent(TRestEvent* evInput) override; void LoadConfig(std::string cfgFilename, std::string name = ""); @@ -51,19 +48,21 @@ class TRestAxionOpticsResponseProcess : public TRestEventProcess { void PrintMetadata() override { BeginPrintProcess(); + RESTMetadata << "Moving to Z: " << fZPosition << " mm" << RESTendl; + EndPrintProcess(); } /// Returns the name of this process - const char* GetProcessName() const override { return "axionOpticsResponse"; } + const char* GetProcessName() const override { return "axionTransport"; } // Constructor - TRestAxionOpticsResponseProcess(); - TRestAxionOpticsResponseProcess(char* cfgFileName); + TRestAxionTransportProcess(); + TRestAxionTransportProcess(char* cfgFileName); // Destructor - ~TRestAxionOpticsResponseProcess(); + ~TRestAxionTransportProcess(); - ClassDefOverride(TRestAxionOpticsResponseProcess, 1); + ClassDefOverride(TRestAxionTransportProcess, 1); }; #endif diff --git a/inc/mpreal.h b/inc/mpreal.h deleted file mode 100644 index 24e286be..00000000 --- a/inc/mpreal.h +++ /dev/null @@ -1,3105 +0,0 @@ -/* - MPFR C++: Multi-precision floating point number class for C++. - Based on MPFR library: http://mpfr.org - - Project homepage: http://www.holoborodko.com/pavel/mpfr - Contact e-mail: pavel@holoborodko.com - - Copyright (c) 2008-2015 Pavel Holoborodko - - Contributors: - Dmitriy Gubanov, Konstantin Holoborodko, Brian Gladman, - Helmut Jarausch, Fokko Beekhof, Ulrich Mutze, Heinz van Saanen, - Pere Constans, Peter van Hoof, Gael Guennebaud, Tsai Chia Cheng, - Alexei Zubanov, Jauhien Piatlicki, Victor Berger, John Westwood, - Petr Aleksandrov, Orion Poplawski, Charles Karney, Arash Partow, - Rodney James, Jorge Leitao. - - Licensing: - (A) MPFR C++ is under GNU General Public License ("GPL"). - - (B) Non-free licenses may also be purchased from the author, for users who - do not want their programs protected by the GPL. - - The non-free licenses are for users that wish to use MPFR C++ in - their products but are unwilling to release their software - under the GPL (which would require them to release source code - and allow free redistribution). - - Such users can purchase an unlimited-use license from the author. - Contact us for more details. - - GNU General Public License ("GPL") copyright permissions statement: - ************************************************************************** - This program 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. - - This program 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 received a copy of the GNU General Public License - along with this program. If not, see . -*/ - -#ifndef __MPREAL_H__ -#define __MPREAL_H__ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// Options -#define MPREAL_HAVE_MSVC_DEBUGVIEW // Enable Debugger Visualizer for "Debug" builds in MSVC. -#define MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS // Enable extended std::numeric_limits - // specialization. -// Meaning that "digits", "round_style" and similar members are defined as functions, not constants. -// See std::numeric_limits at the end of the file for more information. - -// Library version -#define MPREAL_VERSION_MAJOR 3 -#define MPREAL_VERSION_MINOR 6 -#define MPREAL_VERSION_PATCHLEVEL 2 -#define MPREAL_VERSION_STRING "3.6.2" - -// Detect compiler using signatures from http://predef.sourceforge.net/ -#if defined(__GNUC__) && defined(__INTEL_COMPILER) -#define IsInf(x) isinf(x) // Intel ICC compiler on Linux - -#elif defined(_MSC_VER) // Microsoft Visual C++ -#define IsInf(x) (!_finite(x)) - -#else -#define IsInf(x) std::isinf(x) // GNU C/C++ (and/or other compilers), just hope for C99 conformance -#endif - -// A Clang feature extension to determine compiler features. -#ifndef __has_feature -#define __has_feature(x) 0 -#endif - -// Detect support for r-value references (move semantic). Borrowed from Eigen. -#if (__has_feature(cxx_rvalue_references) || defined(__GXX_EXPERIMENTAL_CXX0X__) || \ - __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1600)) - -#define MPREAL_HAVE_MOVE_SUPPORT - -// Use fields in mpfr_t structure to check if it was initialized / set dummy initialization -#define mpfr_is_initialized(x) (0 != (x)->_mpfr_d) -#define mpfr_set_uninitialized(x) ((x)->_mpfr_d = 0) -#endif - -// Detect support for explicit converters. -#if (__has_feature(cxx_explicit_conversions) || \ - (defined(__GXX_EXPERIMENTAL_CXX0X__) && __GNUC_MINOR >= 5) || __cplusplus >= 201103L || \ - (defined(_MSC_VER) && _MSC_VER >= 1800)) - -#define MPREAL_HAVE_EXPLICIT_CONVERTERS -#endif - -#define MPFR_USE_INTMAX_T // Enable 64-bit integer types - should be defined before mpfr.h - -#if defined(MPREAL_HAVE_MSVC_DEBUGVIEW) && defined(_MSC_VER) && defined(_DEBUG) -#define MPREAL_MSVC_DEBUGVIEW_CODE DebugView = toString(); -#define MPREAL_MSVC_DEBUGVIEW_DATA std::string DebugView; -#else -#define MPREAL_MSVC_DEBUGVIEW_CODE -#define MPREAL_MSVC_DEBUGVIEW_DATA -#endif - -#include - -#if (MPFR_VERSION < MPFR_VERSION_NUM(3, 0, 0)) -#include // Needed for random() -#endif - -// Less important options -#define MPREAL_DOUBLE_BITS_OVERFLOW -1 // Triggers overflow exception during conversion to double if mpreal - // cannot fit in MPREAL_DOUBLE_BITS_OVERFLOW bits - // = -1 disables overflow checks (default) - -// Fast replacement for mpfr_set_zero(x, +1): -// (a) uses low-level data members, might not be compatible with new versions of MPFR -// (b) sign is not set, add (x)->_mpfr_sign = 1; -#define mpfr_set_zero_fast(x) ((x)->_mpfr_exp = __MPFR_EXP_ZERO) - -#if defined(__GNUC__) -#define MPREAL_PERMISSIVE_EXPR __extension__ -#else -#define MPREAL_PERMISSIVE_EXPR -#endif - -namespace mpfr { - -class mpreal { - private: - mpfr_t mp; - - public: - // Get default rounding mode & precision - inline static mp_rnd_t get_default_rnd() { return (mp_rnd_t)(mpfr_get_default_rounding_mode()); } - inline static mp_prec_t get_default_prec() { return mpfr_get_default_prec(); } - - // Constructors && type conversions - mpreal(); - mpreal(const mpreal& u); - mpreal(const mpf_t u); - mpreal(const mpz_t u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const mpq_t u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const double u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const long double u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const unsigned long long int u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const long long int u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const unsigned long int u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const unsigned int u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const long int u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const int u, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t mode = mpreal::get_default_rnd()); - - // Construct mpreal from mpfr_t structure. - // shared = true allows to avoid deep copy, so that mpreal and 'u' share the same data & pointers. - mpreal(const mpfr_t u, bool shared = false); - - mpreal(const char* s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, - mp_rnd_t mode = mpreal::get_default_rnd()); - mpreal(const std::string& s, mp_prec_t prec = mpreal::get_default_prec(), int base = 10, - mp_rnd_t mode = mpreal::get_default_rnd()); - - ~mpreal(); - -#ifdef MPREAL_HAVE_MOVE_SUPPORT - mpreal& operator=(mpreal&& v); - mpreal(mpreal&& u); -#endif - - // Operations - // = - // +, -, *, /, ++, --, <<, >> - // *=, +=, -=, /=, - // <, >, ==, <=, >= - - // = - mpreal& operator=(const mpreal& v); - mpreal& operator=(const mpf_t v); - mpreal& operator=(const mpz_t v); - mpreal& operator=(const mpq_t v); - mpreal& operator=(const long double v); - mpreal& operator=(const double v); - mpreal& operator=(const unsigned long int v); - mpreal& operator=(const unsigned long long int v); - mpreal& operator=(const long long int v); - mpreal& operator=(const unsigned int v); - mpreal& operator=(const long int v); - mpreal& operator=(const int v); - mpreal& operator=(const char* s); - mpreal& operator=(const std::string& s); - template - mpreal& operator=(const std::complex& z); - - // + - mpreal& operator+=(const mpreal& v); - mpreal& operator+=(const mpf_t v); - mpreal& operator+=(const mpz_t v); - mpreal& operator+=(const mpq_t v); - mpreal& operator+=(const long double u); - mpreal& operator+=(const double u); - mpreal& operator+=(const unsigned long int u); - mpreal& operator+=(const unsigned int u); - mpreal& operator+=(const long int u); - mpreal& operator+=(const int u); - - mpreal& operator+=(const long long int u); - mpreal& operator+=(const unsigned long long int u); - mpreal& operator-=(const long long int u); - mpreal& operator-=(const unsigned long long int u); - mpreal& operator*=(const long long int u); - mpreal& operator*=(const unsigned long long int u); - mpreal& operator/=(const long long int u); - mpreal& operator/=(const unsigned long long int u); - - const mpreal operator+() const; - mpreal& operator++(); - const mpreal operator++(int); - - // - - mpreal& operator-=(const mpreal& v); - mpreal& operator-=(const mpz_t v); - mpreal& operator-=(const mpq_t v); - mpreal& operator-=(const long double u); - mpreal& operator-=(const double u); - mpreal& operator-=(const unsigned long int u); - mpreal& operator-=(const unsigned int u); - mpreal& operator-=(const long int u); - mpreal& operator-=(const int u); - const mpreal operator-() const; - friend const mpreal operator-(const unsigned long int b, const mpreal& a); - friend const mpreal operator-(const unsigned int b, const mpreal& a); - friend const mpreal operator-(const long int b, const mpreal& a); - friend const mpreal operator-(const int b, const mpreal& a); - friend const mpreal operator-(const double b, const mpreal& a); - mpreal& operator--(); - const mpreal operator--(int); - - // * - mpreal& operator*=(const mpreal& v); - mpreal& operator*=(const mpz_t v); - mpreal& operator*=(const mpq_t v); - mpreal& operator*=(const long double v); - mpreal& operator*=(const double v); - mpreal& operator*=(const unsigned long int v); - mpreal& operator*=(const unsigned int v); - mpreal& operator*=(const long int v); - mpreal& operator*=(const int v); - - // / - mpreal& operator/=(const mpreal& v); - mpreal& operator/=(const mpz_t v); - mpreal& operator/=(const mpq_t v); - mpreal& operator/=(const long double v); - mpreal& operator/=(const double v); - mpreal& operator/=(const unsigned long int v); - mpreal& operator/=(const unsigned int v); - mpreal& operator/=(const long int v); - mpreal& operator/=(const int v); - friend const mpreal operator/(const unsigned long int b, const mpreal& a); - friend const mpreal operator/(const unsigned int b, const mpreal& a); - friend const mpreal operator/(const long int b, const mpreal& a); - friend const mpreal operator/(const int b, const mpreal& a); - friend const mpreal operator/(const double b, const mpreal& a); - - //<<= Fast Multiplication by 2^u - mpreal& operator<<=(const unsigned long int u); - mpreal& operator<<=(const unsigned int u); - mpreal& operator<<=(const long int u); - mpreal& operator<<=(const int u); - - //>>= Fast Division by 2^u - mpreal& operator>>=(const unsigned long int u); - mpreal& operator>>=(const unsigned int u); - mpreal& operator>>=(const long int u); - mpreal& operator>>=(const int u); - - // Type Conversion operators - bool toBool() const; - long toLong(mp_rnd_t mode = GMP_RNDZ) const; - unsigned long toULong(mp_rnd_t mode = GMP_RNDZ) const; - long long toLLong(mp_rnd_t mode = GMP_RNDZ) const; - unsigned long long toULLong(mp_rnd_t mode = GMP_RNDZ) const; - float toFloat(mp_rnd_t mode = GMP_RNDN) const; - double toDouble(mp_rnd_t mode = GMP_RNDN) const; - long double toLDouble(mp_rnd_t mode = GMP_RNDN) const; - -#if defined(MPREAL_HAVE_EXPLICIT_CONVERTERS) - explicit operator bool() const { return toBool(); } - explicit operator int() const { return int(toLong()); } - explicit operator long() const { return toLong(); } - explicit operator long long() const { return toLLong(); } - explicit operator unsigned() const { return unsigned(toULong()); } - explicit operator unsigned long() const { return toULong(); } - explicit operator unsigned long long() const { return toULLong(); } - explicit operator float() const { return toFloat(); } - explicit operator double() const { return toDouble(); } - explicit operator long double() const { return toLDouble(); } -#endif - - // Get raw pointers so that mpreal can be directly used in raw mpfr_* functions - ::mpfr_ptr mpfr_ptr(); - ::mpfr_srcptr mpfr_ptr() const; - ::mpfr_srcptr mpfr_srcptr() const; - - // Convert mpreal to string with n significant digits in base b - // n = -1 -> convert with the maximum available digits - std::string toString(int n = -1, int b = 10, mp_rnd_t mode = mpreal::get_default_rnd()) const; - -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - std::string toString(const std::string& format) const; -#endif - - std::ostream& output(std::ostream& os) const; - - // Math Functions - friend const mpreal sqr(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal sqrt(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal sqrt(const unsigned long int v, mp_rnd_t rnd_mode); - friend const mpreal cbrt(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal root(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode); - friend const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode); - friend const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode); - friend const mpreal pow(const mpreal& a, const unsigned long int b, mp_rnd_t rnd_mode); - friend const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode); - friend const mpreal pow(const unsigned long int a, const mpreal& b, mp_rnd_t rnd_mode); - friend const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode); - friend const mpreal fabs(const mpreal& v, mp_rnd_t rnd_mode); - - friend const mpreal abs(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode); - friend inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode); - friend inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode); - friend inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode); - friend inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode); - friend int cmpabs(const mpreal& a, const mpreal& b); - - friend const mpreal log(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal log2(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal logb(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal log10(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal exp(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal exp2(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal exp10(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal log1p(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal expm1(const mpreal& v, mp_rnd_t rnd_mode); - - friend const mpreal cos(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal sin(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal tan(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal sec(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal csc(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal cot(const mpreal& v, mp_rnd_t rnd_mode); - friend int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode); - - friend const mpreal acos(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal asin(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal atan(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal atan2(const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode); - friend const mpreal acot(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal asec(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal acsc(const mpreal& v, mp_rnd_t rnd_mode); - - friend const mpreal cosh(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal sinh(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal tanh(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal sech(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal csch(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal coth(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal acosh(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal asinh(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal atanh(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal acoth(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal asech(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal acsch(const mpreal& v, mp_rnd_t rnd_mode); - - friend const mpreal hypot(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); - - friend const mpreal fac_ui(unsigned long int v, mp_prec_t prec, mp_rnd_t rnd_mode); - friend const mpreal eint(const mpreal& v, mp_rnd_t rnd_mode); - - friend const mpreal gamma(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal tgamma(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal lngamma(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal lgamma(const mpreal& v, int* signp, mp_rnd_t rnd_mode); - friend const mpreal zeta(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal erf(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal erfc(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal besselj0(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal besselj1(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal besseljn(long n, const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal bessely0(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal bessely1(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal besselyn(long n, const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal fma(const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode); - friend const mpreal fms(const mpreal& v1, const mpreal& v2, const mpreal& v3, mp_rnd_t rnd_mode); - friend const mpreal agm(const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode); - friend const mpreal sum(const mpreal tab[], const unsigned long int n, int& status, mp_rnd_t rnd_mode); - friend int sgn(const mpreal& v); // returns -1 or +1 - -// MPFR 2.4.0 Specifics -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - friend int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal li2(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal fmod(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); - friend const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode); - - // MATLAB's semantic equivalents - friend const mpreal rem(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Remainder after division - friend const mpreal mod(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); // Modulus after division -#endif - -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3, 0, 0)) - friend const mpreal digamma(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal ai(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal urandom( - gmp_randstate_t& state, - mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear -#endif - -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3, 1, 0)) - friend const mpreal grandom( - gmp_randstate_t& state, - mp_rnd_t rnd_mode); // use gmp_randinit_default() to init state, gmp_randclear() to clear - friend const mpreal grandom(unsigned int seed); -#endif - - // Uniformly distributed random number generation in [0,1] using - // Mersenne-Twister algorithm by default. - // Use parameter to setup seed, e.g.: random((unsigned)time(NULL)) - // Check urandom() for more precise control. - friend const mpreal random(unsigned int seed); - - // Exponent and mantissa manipulation - friend const mpreal frexp(const mpreal& v, mp_exp_t* exp); - friend const mpreal ldexp(const mpreal& v, mp_exp_t exp); - friend const mpreal scalbn(const mpreal& v, mp_exp_t exp); - - // Splits mpreal value into fractional and integer parts. - // Returns fractional part and stores integer part in n. - friend const mpreal modf(const mpreal& v, mpreal& n); - - // Constants - // don't forget to call mpfr_free_cache() for every thread where you are using const-functions - friend const mpreal const_log2(mp_prec_t prec, mp_rnd_t rnd_mode); - friend const mpreal const_pi(mp_prec_t prec, mp_rnd_t rnd_mode); - friend const mpreal const_euler(mp_prec_t prec, mp_rnd_t rnd_mode); - friend const mpreal const_catalan(mp_prec_t prec, mp_rnd_t rnd_mode); - - // returns +inf iff sign>=0 otherwise -inf - friend const mpreal const_infinity(int sign, mp_prec_t prec); - - // Output/ Input - friend std::ostream& operator<<(std::ostream& os, const mpreal& v); - friend std::istream& operator>>(std::istream& is, mpreal& v); - - // Integer Related Functions - friend const mpreal rint(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal ceil(const mpreal& v); - friend const mpreal floor(const mpreal& v); - friend const mpreal round(const mpreal& v); - friend const mpreal trunc(const mpreal& v); - friend const mpreal rint_ceil(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal rint_floor(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal rint_round(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal rint_trunc(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal frac(const mpreal& v, mp_rnd_t rnd_mode); - friend const mpreal remainder(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); - friend const mpreal remquo(long* q, const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); - - // Miscellaneous Functions - friend const mpreal nexttoward(const mpreal& x, const mpreal& y); - friend const mpreal nextabove(const mpreal& x); - friend const mpreal nextbelow(const mpreal& x); - - // use gmp_randinit_default() to init state, gmp_randclear() to clear - friend const mpreal urandomb(gmp_randstate_t& state); - -// MPFR < 2.4.2 Specifics -#if (MPFR_VERSION <= MPFR_VERSION_NUM(2, 4, 2)) - friend const mpreal random2(mp_size_t size, mp_exp_t exp); -#endif - - // Instance Checkers - friend bool isnan(const mpreal& v); - friend bool isinf(const mpreal& v); - friend bool isfinite(const mpreal& v); - - friend bool isnum(const mpreal& v); - friend bool iszero(const mpreal& v); - friend bool isint(const mpreal& v); - -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3, 0, 0)) - friend bool isregular(const mpreal& v); -#endif - - // Set/Get instance properties - inline mp_prec_t get_prec() const; - inline void set_prec(mp_prec_t prec, - mp_rnd_t rnd_mode = get_default_rnd()); // Change precision with rounding mode - - // Aliases for get_prec(), set_prec() - needed for compatibility with std::complex interface - inline mpreal& setPrecision(int Precision, mp_rnd_t RoundingMode = get_default_rnd()); - inline int getPrecision() const; - - // Set mpreal to +/- inf, NaN, +/-0 - mpreal& setInf(int Sign = +1); - mpreal& setNan(); - mpreal& setZero(int Sign = +1); - mpreal& setSign(int Sign, mp_rnd_t RoundingMode = get_default_rnd()); - - // Exponent - mp_exp_t get_exp(); - int set_exp(mp_exp_t e); - int check_range(int t, mp_rnd_t rnd_mode = get_default_rnd()); - int subnormalize(int t, mp_rnd_t rnd_mode = get_default_rnd()); - - // Inexact conversion from float - inline bool fits_in_bits(double x, int n); - - // Set/Get global properties - static void set_default_prec(mp_prec_t prec); - static void set_default_rnd(mp_rnd_t rnd_mode); - - static mp_exp_t get_emin(void); - static mp_exp_t get_emax(void); - static mp_exp_t get_emin_min(void); - static mp_exp_t get_emin_max(void); - static mp_exp_t get_emax_min(void); - static mp_exp_t get_emax_max(void); - static int set_emin(mp_exp_t exp); - static int set_emax(mp_exp_t exp); - - // Efficient swapping of two mpreal values - needed for std algorithms - friend void swap(mpreal& x, mpreal& y); - - friend const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); - friend const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode); - - private: - // Human friendly Debug Preview in Visual Studio. - // Put one of these lines: - // - // mpfr::mpreal= ; Show value only - // mpfr::mpreal=, bits ; Show value & precision - // - // at the beginning of - // [Visual Studio Installation Folder]\Common7\Packages\Debugger\autoexp.dat - MPREAL_MSVC_DEBUGVIEW_DATA - - // "Smart" resources deallocation. Checks if instance initialized before deletion. - void clear(::mpfr_ptr); -}; - -////////////////////////////////////////////////////////////////////////// -// Exceptions -class conversion_overflow : public std::exception { - public: - std::string why() { return "inexact conversion from floating point"; } -}; - -////////////////////////////////////////////////////////////////////////// -// Constructors & converters -// Default constructor: creates mp number and initializes it to 0. -inline mpreal::mpreal() { - mpfr_init2(mpfr_ptr(), mpreal::get_default_prec()); - mpfr_set_zero_fast(mpfr_ptr()); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const mpreal& u) { - mpfr_init2(mpfr_ptr(), mpfr_get_prec(u.mpfr_srcptr())); - mpfr_set(mpfr_ptr(), u.mpfr_srcptr(), mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -#ifdef MPREAL_HAVE_MOVE_SUPPORT -inline mpreal::mpreal(mpreal&& other) { - mpfr_set_uninitialized(mpfr_ptr()); // make sure "other" holds no pointer to actual data - mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal& mpreal::operator=(mpreal&& other) { - mpfr_swap(mpfr_ptr(), other.mpfr_ptr()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} -#endif - -inline mpreal::mpreal(const mpfr_t u, bool shared) { - if (shared) { - std::memcpy(mpfr_ptr(), u, sizeof(mpfr_t)); - } else { - mpfr_init2(mpfr_ptr(), mpfr_get_prec(u)); - mpfr_set(mpfr_ptr(), u, mpreal::get_default_rnd()); - } - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const mpf_t u) { - mpfr_init2(mpfr_ptr(), - (mp_prec_t)mpf_get_prec(u)); // (gmp: mp_bitcnt_t) unsigned long -> long (mpfr: mp_prec_t) - mpfr_set_f(mpfr_ptr(), u, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const mpz_t u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_z(mpfr_ptr(), u, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const mpq_t u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_q(mpfr_ptr(), u, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const double u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - -#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1) - if (fits_in_bits(u, MPREAL_DOUBLE_BITS_OVERFLOW)) { - mpfr_set_d(mpfr_ptr(), u, mode); - } else - throw conversion_overflow(); -#else - mpfr_set_d(mpfr_ptr(), u, mode); -#endif - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const long double u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_ld(mpfr_ptr(), u, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const unsigned long long int u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_uj(mpfr_ptr(), u, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const long long int u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_sj(mpfr_ptr(), u, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const unsigned long int u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_ui(mpfr_ptr(), u, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const unsigned int u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_ui(mpfr_ptr(), u, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const long int u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_si(mpfr_ptr(), u, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const int u, mp_prec_t prec, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_si(mpfr_ptr(), u, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_str(mpfr_ptr(), s, base, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mpreal::mpreal(const std::string& s, mp_prec_t prec, int base, mp_rnd_t mode) { - mpfr_init2(mpfr_ptr(), prec); - mpfr_set_str(mpfr_ptr(), s.c_str(), base, mode); - - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline void mpreal::clear(::mpfr_ptr x) { -#ifdef MPREAL_HAVE_MOVE_SUPPORT - if (mpfr_is_initialized(x)) -#endif - mpfr_clear(x); -} - -inline mpreal::~mpreal() { clear(mpfr_ptr()); } - -// internal namespace needed for template magic -namespace internal { - -// Use SFINAE to restrict arithmetic operations instantiation only for numeric types -// This is needed for smooth integration with libraries based on expression templates, like Eigen. -// TODO: Do the same for boolean operators. -template -struct result_type {}; - -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -template <> -struct result_type { - typedef mpreal type; -}; -} // namespace internal - -// + Addition -template -inline const typename internal::result_type::type operator+(const mpreal& lhs, const Rhs& rhs) { - return mpreal(lhs) += rhs; -} - -template -inline const typename internal::result_type::type operator+(const Lhs& lhs, const mpreal& rhs) { - return mpreal(rhs) += lhs; -} - -// - Subtraction -template -inline const typename internal::result_type::type operator-(const mpreal& lhs, const Rhs& rhs) { - return mpreal(lhs) -= rhs; -} - -template -inline const typename internal::result_type::type operator-(const Lhs& lhs, const mpreal& rhs) { - return mpreal(lhs) -= rhs; -} - -// * Multiplication -template -inline const typename internal::result_type::type operator*(const mpreal& lhs, const Rhs& rhs) { - return mpreal(lhs) *= rhs; -} - -template -inline const typename internal::result_type::type operator*(const Lhs& lhs, const mpreal& rhs) { - return mpreal(rhs) *= lhs; -} - -// / Division -template -inline const typename internal::result_type::type operator/(const mpreal& lhs, const Rhs& rhs) { - return mpreal(lhs) /= rhs; -} - -template -inline const typename internal::result_type::type operator/(const Lhs& lhs, const mpreal& rhs) { - return mpreal(lhs) /= rhs; -} - -////////////////////////////////////////////////////////////////////////// -// sqrt -const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal sqrt(const long int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal sqrt(const int v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal sqrt(const long double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal sqrt(const double v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -// abs -inline const mpreal abs(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()); - -////////////////////////////////////////////////////////////////////////// -// pow -const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -const mpreal pow(const unsigned long int a, const unsigned int b, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const unsigned long int a, const long double b, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -const mpreal pow(const unsigned int a, const unsigned long int b, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long double a, const unsigned long int b, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -inline const mpreal div_2ui(const mpreal& v, unsigned long int k, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()); -inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode = mpreal::get_default_rnd()); - -////////////////////////////////////////////////////////////////////////// -// Estimate machine epsilon for the given precision -// Returns smallest eps such that 1.0 + eps != 1.0 -inline mpreal machine_epsilon(mp_prec_t prec = mpreal::get_default_prec()); - -// Returns smallest eps such that x + eps != x (relative machine epsilon) -inline mpreal machine_epsilon(const mpreal& x); - -// Gives max & min values for the required precision, -// minval is 'safe' meaning 1 / minval does not overflow -// maxval is 'safe' meaning 1 / maxval does not underflow -inline mpreal minval(mp_prec_t prec = mpreal::get_default_prec()); -inline mpreal maxval(mp_prec_t prec = mpreal::get_default_prec()); - -// 'Dirty' equality check 1: |a-b| < min{|a|,|b|} * eps -inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps); - -// 'Dirty' equality check 2: |a-b| < min{|a|,|b|} * eps( min{|a|,|b|} ) -inline bool isEqualFuzzy(const mpreal& a, const mpreal& b); - -// 'Bitwise' equality check -// maxUlps - a and b can be apart by maxUlps binary numbers. -inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps); - -////////////////////////////////////////////////////////////////////////// -// Convert precision in 'bits' to decimal digits and vice versa. -// bits = ceil(digits*log[2](10)) -// digits = floor(bits*log[10](2)) - -inline mp_prec_t digits2bits(int d); -inline int bits2digits(mp_prec_t b); - -////////////////////////////////////////////////////////////////////////// -// min, max -const mpreal(max)(const mpreal& x, const mpreal& y); -const mpreal(min)(const mpreal& x, const mpreal& y); - -////////////////////////////////////////////////////////////////////////// -// Implementation -////////////////////////////////////////////////////////////////////////// - -////////////////////////////////////////////////////////////////////////// -// Operators - Assignment -inline mpreal& mpreal::operator=(const mpreal& v) { - if (this != &v) { - mp_prec_t tp = mpfr_get_prec(mpfr_srcptr()); - mp_prec_t vp = mpfr_get_prec(v.mpfr_srcptr()); - - if (tp != vp) { - clear(mpfr_ptr()); - mpfr_init2(mpfr_ptr(), vp); - } - - mpfr_set(mpfr_ptr(), v.mpfr_srcptr(), mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - } - return *this; -} - -inline mpreal& mpreal::operator=(const mpf_t v) { - mpfr_set_f(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const mpz_t v) { - mpfr_set_z(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const mpq_t v) { - mpfr_set_q(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const long double v) { - mpfr_set_ld(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const double v) { -#if (MPREAL_DOUBLE_BITS_OVERFLOW > -1) - if (fits_in_bits(v, MPREAL_DOUBLE_BITS_OVERFLOW)) { - mpfr_set_d(mpfr_ptr(), v, mpreal::get_default_rnd()); - } else - throw conversion_overflow(); -#else - mpfr_set_d(mpfr_ptr(), v, mpreal::get_default_rnd()); -#endif - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const unsigned long int v) { - mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const unsigned int v) { - mpfr_set_ui(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const unsigned long long int v) { - mpfr_set_uj(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const long long int v) { - mpfr_set_sj(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const long int v) { - mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const int v) { - mpfr_set_si(mpfr_ptr(), v, mpreal::get_default_rnd()); - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator=(const char* s) { - // Use other converters for more precise control on base & precision & rounding: - // - // mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) - // mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode) - // - // Here we assume base = 10 and we use precision of target variable. - - mpfr_t t; - - mpfr_init2(t, mpfr_get_prec(mpfr_srcptr())); - - if (0 == mpfr_set_str(t, s, 10, mpreal::get_default_rnd())) { - mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - } - - clear(t); - return *this; -} - -inline mpreal& mpreal::operator=(const std::string& s) { - // Use other converters for more precise control on base & precision & rounding: - // - // mpreal(const char* s, mp_prec_t prec, int base, mp_rnd_t mode) - // mpreal(const std::string& s,mp_prec_t prec, int base, mp_rnd_t mode) - // - // Here we assume base = 10 and we use precision of target variable. - - mpfr_t t; - - mpfr_init2(t, mpfr_get_prec(mpfr_srcptr())); - - if (0 == mpfr_set_str(t, s.c_str(), 10, mpreal::get_default_rnd())) { - mpfr_set(mpfr_ptr(), t, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - } - - clear(t); - return *this; -} - -template -inline mpreal& mpreal::operator=(const std::complex& z) { - return *this = z.real(); -} - -////////////////////////////////////////////////////////////////////////// -// + Addition -inline mpreal& mpreal::operator+=(const mpreal& v) { - mpfr_add(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const mpf_t u) { - *this += mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const mpz_t u) { - mpfr_add_z(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const mpq_t u) { - mpfr_add_q(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const long double u) { - *this += mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const double u) { -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - mpfr_add_d(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); -#else - *this += mpreal(u); -#endif - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const unsigned long int u) { - mpfr_add_ui(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const unsigned int u) { - mpfr_add_ui(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const long int u) { - mpfr_add_si(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const int u) { - mpfr_add_si(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator+=(const long long int u) { - *this += mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} -inline mpreal& mpreal::operator+=(const unsigned long long int u) { - *this += mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} -inline mpreal& mpreal::operator-=(const long long int u) { - *this -= mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} -inline mpreal& mpreal::operator-=(const unsigned long long int u) { - *this -= mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} -inline mpreal& mpreal::operator*=(const long long int u) { - *this *= mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} -inline mpreal& mpreal::operator*=(const unsigned long long int u) { - *this *= mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} -inline mpreal& mpreal::operator/=(const long long int u) { - *this /= mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} -inline mpreal& mpreal::operator/=(const unsigned long long int u) { - *this /= mpreal(u); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline const mpreal mpreal::operator+() const { return mpreal(*this); } - -inline const mpreal operator+(const mpreal& a, const mpreal& b) { - mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); - mpfr_add(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); - return c; -} - -inline mpreal& mpreal::operator++() { return *this += 1; } - -inline const mpreal mpreal::operator++(int) { - mpreal x(*this); - *this += 1; - return x; -} - -inline mpreal& mpreal::operator--() { return *this -= 1; } - -inline const mpreal mpreal::operator--(int) { - mpreal x(*this); - *this -= 1; - return x; -} - -////////////////////////////////////////////////////////////////////////// -// - Subtraction -inline mpreal& mpreal::operator-=(const mpreal& v) { - mpfr_sub(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator-=(const mpz_t v) { - mpfr_sub_z(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator-=(const mpq_t v) { - mpfr_sub_q(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator-=(const long double v) { - *this -= mpreal(v); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator-=(const double v) { -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - mpfr_sub_d(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); -#else - *this -= mpreal(v); -#endif - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator-=(const unsigned long int v) { - mpfr_sub_ui(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator-=(const unsigned int v) { - mpfr_sub_ui(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator-=(const long int v) { - mpfr_sub_si(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator-=(const int v) { - mpfr_sub_si(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline const mpreal mpreal::operator-() const { - mpreal u(*this); - mpfr_neg(u.mpfr_ptr(), u.mpfr_srcptr(), mpreal::get_default_rnd()); - return u; -} - -inline const mpreal operator-(const mpreal& a, const mpreal& b) { - mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); - mpfr_sub(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); - return c; -} - -inline const mpreal operator-(const double b, const mpreal& a) { -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); - mpfr_d_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -#else - mpreal x(b, mpfr_get_prec(a.mpfr_ptr())); - x -= a; - return x; -#endif -} - -inline const mpreal operator-(const unsigned long int b, const mpreal& a) { - mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); - mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -} - -inline const mpreal operator-(const unsigned int b, const mpreal& a) { - mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); - mpfr_ui_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -} - -inline const mpreal operator-(const long int b, const mpreal& a) { - mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); - mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -} - -inline const mpreal operator-(const int b, const mpreal& a) { - mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); - mpfr_si_sub(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -} - -////////////////////////////////////////////////////////////////////////// -// * Multiplication -inline mpreal& mpreal::operator*=(const mpreal& v) { - mpfr_mul(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator*=(const mpz_t v) { - mpfr_mul_z(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator*=(const mpq_t v) { - mpfr_mul_q(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator*=(const long double v) { - *this *= mpreal(v); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator*=(const double v) { -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - mpfr_mul_d(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); -#else - *this *= mpreal(v); -#endif - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator*=(const unsigned long int v) { - mpfr_mul_ui(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator*=(const unsigned int v) { - mpfr_mul_ui(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator*=(const long int v) { - mpfr_mul_si(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator*=(const int v) { - mpfr_mul_si(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline const mpreal operator*(const mpreal& a, const mpreal& b) { - mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_ptr()), mpfr_get_prec(b.mpfr_ptr()))); - mpfr_mul(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); - return c; -} - -////////////////////////////////////////////////////////////////////////// -// / Division -inline mpreal& mpreal::operator/=(const mpreal& v) { - mpfr_div(mpfr_ptr(), mpfr_srcptr(), v.mpfr_srcptr(), mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator/=(const mpz_t v) { - mpfr_div_z(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator/=(const mpq_t v) { - mpfr_div_q(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator/=(const long double v) { - *this /= mpreal(v); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator/=(const double v) { -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - mpfr_div_d(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); -#else - *this /= mpreal(v); -#endif - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator/=(const unsigned long int v) { - mpfr_div_ui(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator/=(const unsigned int v) { - mpfr_div_ui(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator/=(const long int v) { - mpfr_div_si(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator/=(const int v) { - mpfr_div_si(mpfr_ptr(), mpfr_srcptr(), v, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline const mpreal operator/(const mpreal& a, const mpreal& b) { - mpreal c(0, (std::max)(mpfr_get_prec(a.mpfr_srcptr()), mpfr_get_prec(b.mpfr_srcptr()))); - mpfr_div(c.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), mpreal::get_default_rnd()); - return c; -} - -inline const mpreal operator/(const unsigned long int b, const mpreal& a) { - mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); - mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -} - -inline const mpreal operator/(const unsigned int b, const mpreal& a) { - mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); - mpfr_ui_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -} - -inline const mpreal operator/(const long int b, const mpreal& a) { - mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); - mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -} - -inline const mpreal operator/(const int b, const mpreal& a) { - mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); - mpfr_si_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -} - -inline const mpreal operator/(const double b, const mpreal& a) { -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - mpreal x(0, mpfr_get_prec(a.mpfr_srcptr())); - mpfr_d_div(x.mpfr_ptr(), b, a.mpfr_srcptr(), mpreal::get_default_rnd()); - return x; -#else - mpreal x(0, mpfr_get_prec(a.mpfr_ptr())); - x /= a; - return x; -#endif -} - -////////////////////////////////////////////////////////////////////////// -// Shifts operators - Multiplication/Division by power of 2 -inline mpreal& mpreal::operator<<=(const unsigned long int u) { - mpfr_mul_2ui(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator<<=(const unsigned int u) { - mpfr_mul_2ui(mpfr_ptr(), mpfr_srcptr(), static_cast(u), mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator<<=(const long int u) { - mpfr_mul_2si(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator<<=(const int u) { - mpfr_mul_2si(mpfr_ptr(), mpfr_srcptr(), static_cast(u), mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator>>=(const unsigned long int u) { - mpfr_div_2ui(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator>>=(const unsigned int u) { - mpfr_div_2ui(mpfr_ptr(), mpfr_srcptr(), static_cast(u), mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator>>=(const long int u) { - mpfr_div_2si(mpfr_ptr(), mpfr_srcptr(), u, mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::operator>>=(const int u) { - mpfr_div_2si(mpfr_ptr(), mpfr_srcptr(), static_cast(u), mpreal::get_default_rnd()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline const mpreal operator<<(const mpreal& v, const unsigned long int k) { return mul_2ui(v, k); } - -inline const mpreal operator<<(const mpreal& v, const unsigned int k) { - return mul_2ui(v, static_cast(k)); -} - -inline const mpreal operator<<(const mpreal& v, const long int k) { return mul_2si(v, k); } - -inline const mpreal operator<<(const mpreal& v, const int k) { return mul_2si(v, static_cast(k)); } - -inline const mpreal operator>>(const mpreal& v, const unsigned long int k) { return div_2ui(v, k); } - -inline const mpreal operator>>(const mpreal& v, const long int k) { return div_2si(v, k); } - -inline const mpreal operator>>(const mpreal& v, const unsigned int k) { - return div_2ui(v, static_cast(k)); -} - -inline const mpreal operator>>(const mpreal& v, const int k) { return div_2si(v, static_cast(k)); } - -// mul_2ui -inline const mpreal mul_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) { - mpreal x(v); - mpfr_mul_2ui(x.mpfr_ptr(), v.mpfr_srcptr(), k, rnd_mode); - return x; -} - -// mul_2si -inline const mpreal mul_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode) { - mpreal x(v); - mpfr_mul_2si(x.mpfr_ptr(), v.mpfr_srcptr(), k, rnd_mode); - return x; -} - -inline const mpreal div_2ui(const mpreal& v, unsigned long int k, mp_rnd_t rnd_mode) { - mpreal x(v); - mpfr_div_2ui(x.mpfr_ptr(), v.mpfr_srcptr(), k, rnd_mode); - return x; -} - -inline const mpreal div_2si(const mpreal& v, long int k, mp_rnd_t rnd_mode) { - mpreal x(v); - mpfr_div_2si(x.mpfr_ptr(), v.mpfr_srcptr(), k, rnd_mode); - return x; -} - -////////////////////////////////////////////////////////////////////////// -// Relational operators - -// WARNING: -// -// Please note that following checks for double-NaN are guaranteed to work only in IEEE math mode: -// -// isnan(b) = (b != b) -// isnan(b) = !(b == b) (we use in code below) -// -// Be cautions if you use compiler options which break strict IEEE compliance (e.g. -ffast-math in GCC). -// Use std::isnan instead (C++11). - -inline bool operator>(const mpreal& a, const mpreal& b) { - return (mpfr_greater_p(a.mpfr_srcptr(), b.mpfr_srcptr()) != 0); -} -inline bool operator>(const mpreal& a, const unsigned long int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) > 0); -} -inline bool operator>(const mpreal& a, const unsigned int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) > 0); -} -inline bool operator>(const mpreal& a, const long int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) > 0); -} -inline bool operator>(const mpreal& a, const int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) > 0); -} -inline bool operator>(const mpreal& a, const long double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(), b) > 0); -} -inline bool operator>(const mpreal& a, const double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_d(a.mpfr_srcptr(), b) > 0); -} - -inline bool operator>=(const mpreal& a, const mpreal& b) { - return (mpfr_greaterequal_p(a.mpfr_srcptr(), b.mpfr_srcptr()) != 0); -} -inline bool operator>=(const mpreal& a, const unsigned long int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) >= 0); -} -inline bool operator>=(const mpreal& a, const unsigned int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) >= 0); -} -inline bool operator>=(const mpreal& a, const long int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) >= 0); -} -inline bool operator>=(const mpreal& a, const int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) >= 0); -} -inline bool operator>=(const mpreal& a, const long double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(), b) >= 0); -} -inline bool operator>=(const mpreal& a, const double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_d(a.mpfr_srcptr(), b) >= 0); -} - -inline bool operator<(const mpreal& a, const mpreal& b) { - return (mpfr_less_p(a.mpfr_srcptr(), b.mpfr_srcptr()) != 0); -} -inline bool operator<(const mpreal& a, const unsigned long int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) < 0); -} -inline bool operator<(const mpreal& a, const unsigned int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) < 0); -} -inline bool operator<(const mpreal& a, const long int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) < 0); -} -inline bool operator<(const mpreal& a, const int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) < 0); -} -inline bool operator<(const mpreal& a, const long double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(), b) < 0); -} -inline bool operator<(const mpreal& a, const double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_d(a.mpfr_srcptr(), b) < 0); -} - -inline bool operator<=(const mpreal& a, const mpreal& b) { - return (mpfr_lessequal_p(a.mpfr_srcptr(), b.mpfr_srcptr()) != 0); -} -inline bool operator<=(const mpreal& a, const unsigned long int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) <= 0); -} -inline bool operator<=(const mpreal& a, const unsigned int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) <= 0); -} -inline bool operator<=(const mpreal& a, const long int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) <= 0); -} -inline bool operator<=(const mpreal& a, const int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) <= 0); -} -inline bool operator<=(const mpreal& a, const long double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(), b) <= 0); -} -inline bool operator<=(const mpreal& a, const double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_d(a.mpfr_srcptr(), b) <= 0); -} - -inline bool operator==(const mpreal& a, const mpreal& b) { - return (mpfr_equal_p(a.mpfr_srcptr(), b.mpfr_srcptr()) != 0); -} -inline bool operator==(const mpreal& a, const unsigned long int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) == 0); -} -inline bool operator==(const mpreal& a, const unsigned int b) { - return !isnan(a) && (mpfr_cmp_ui(a.mpfr_srcptr(), b) == 0); -} -inline bool operator==(const mpreal& a, const long int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) == 0); -} -inline bool operator==(const mpreal& a, const int b) { - return !isnan(a) && (mpfr_cmp_si(a.mpfr_srcptr(), b) == 0); -} -inline bool operator==(const mpreal& a, const long double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_ld(a.mpfr_srcptr(), b) == 0); -} -inline bool operator==(const mpreal& a, const double b) { - return !isnan(a) && (b == b) && (mpfr_cmp_d(a.mpfr_srcptr(), b) == 0); -} - -inline bool operator!=(const mpreal& a, const mpreal& b) { return !(a == b); } -inline bool operator!=(const mpreal& a, const unsigned long int b) { return !(a == b); } -inline bool operator!=(const mpreal& a, const unsigned int b) { return !(a == b); } -inline bool operator!=(const mpreal& a, const long int b) { return !(a == b); } -inline bool operator!=(const mpreal& a, const int b) { return !(a == b); } -inline bool operator!=(const mpreal& a, const long double b) { return !(a == b); } -inline bool operator!=(const mpreal& a, const double b) { return !(a == b); } - -inline bool isnan(const mpreal& op) { return (mpfr_nan_p(op.mpfr_srcptr()) != 0); } -inline bool isinf(const mpreal& op) { return (mpfr_inf_p(op.mpfr_srcptr()) != 0); } -inline bool isfinite(const mpreal& op) { return (mpfr_number_p(op.mpfr_srcptr()) != 0); } -inline bool iszero(const mpreal& op) { return (mpfr_zero_p(op.mpfr_srcptr()) != 0); } -inline bool isint(const mpreal& op) { return (mpfr_integer_p(op.mpfr_srcptr()) != 0); } - -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3, 0, 0)) -inline bool isregular(const mpreal& op) { return (mpfr_regular_p(op.mpfr_srcptr())); } -#endif - -////////////////////////////////////////////////////////////////////////// -// Type Converters -inline bool mpreal::toBool() const { return mpfr_zero_p(mpfr_srcptr()) == 0; } -inline long mpreal::toLong(mp_rnd_t mode) const { return mpfr_get_si(mpfr_srcptr(), mode); } -inline unsigned long mpreal::toULong(mp_rnd_t mode) const { return mpfr_get_ui(mpfr_srcptr(), mode); } -inline float mpreal::toFloat(mp_rnd_t mode) const { return mpfr_get_flt(mpfr_srcptr(), mode); } -inline double mpreal::toDouble(mp_rnd_t mode) const { return mpfr_get_d(mpfr_srcptr(), mode); } -inline long double mpreal::toLDouble(mp_rnd_t mode) const { return mpfr_get_ld(mpfr_srcptr(), mode); } -inline long long mpreal::toLLong(mp_rnd_t mode) const { return mpfr_get_sj(mpfr_srcptr(), mode); } -inline unsigned long long mpreal::toULLong(mp_rnd_t mode) const { return mpfr_get_uj(mpfr_srcptr(), mode); } - -inline ::mpfr_ptr mpreal::mpfr_ptr() { return mp; } -inline ::mpfr_srcptr mpreal::mpfr_ptr() const { return mp; } -inline ::mpfr_srcptr mpreal::mpfr_srcptr() const { return mp; } - -template -inline std::string toString(T t, std::ios_base& (*f)(std::ios_base&)) { - std::ostringstream oss; - oss << f << t; - return oss.str(); -} - -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - -inline std::string mpreal::toString(const std::string& format) const { - char* s = NULL; - std::string out; - - if (!format.empty()) { - if (!(mpfr_asprintf(&s, format.c_str(), mpfr_srcptr()) < 0)) { - out = std::string(s); - - mpfr_free_str(s); - } - } - - return out; -} - -#endif - -inline std::string mpreal::toString(int n, int b, mp_rnd_t mode) const { - // TODO: Add extended format specification (f, e, rounding mode) as it done in output operator - (void)b; - (void)mode; - -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - - std::ostringstream format; - - int digits = (n >= 0) ? n : 1 + bits2digits(mpfr_get_prec(mpfr_srcptr())); - - format << "%." << digits << "RNg"; - - return toString(format.str()); - -#else - - char *s, *ns = NULL; - size_t slen, nslen; - mp_exp_t exp; - std::string out; - - if (mpfr_inf_p(mp)) { - if (mpfr_sgn(mp) > 0) - return "+Inf"; - else - return "-Inf"; - } - - if (mpfr_zero_p(mp)) return "0"; - if (mpfr_nan_p(mp)) return "NaN"; - - s = mpfr_get_str(NULL, &exp, b, 0, mp, mode); - ns = mpfr_get_str(NULL, &exp, b, (std::max)(0, n), mp, mode); - - if (s != NULL && ns != NULL) { - slen = strlen(s); - nslen = strlen(ns); - if (nslen <= slen) { - mpfr_free_str(s); - s = ns; - slen = nslen; - } else { - mpfr_free_str(ns); - } - - // Make human eye-friendly formatting if possible - if (exp > 0 && static_cast(exp) < slen) { - if (s[0] == '-') { - // Remove zeros starting from right end - char* ptr = s + slen - 1; - while (*ptr == '0' && ptr > s + exp) ptr--; - - if (ptr == s + exp) - out = std::string(s, exp + 1); - else - out = std::string(s, exp + 1) + '.' + std::string(s + exp + 1, ptr - (s + exp + 1) + 1); - - // out = string(s,exp+1)+'.'+string(s+exp+1); - } else { - // Remove zeros starting from right end - char* ptr = s + slen - 1; - while (*ptr == '0' && ptr > s + exp - 1) ptr--; - - if (ptr == s + exp - 1) - out = std::string(s, exp); - else - out = std::string(s, exp) + '.' + std::string(s + exp, ptr - (s + exp) + 1); - - // out = string(s,exp)+'.'+string(s+exp); - } - - } else { // exp<0 || exp>slen - if (s[0] == '-') { - // Remove zeros starting from right end - char* ptr = s + slen - 1; - while (*ptr == '0' && ptr > s + 1) ptr--; - - if (ptr == s + 1) - out = std::string(s, 2); - else - out = std::string(s, 2) + '.' + std::string(s + 2, ptr - (s + 2) + 1); - - // out = string(s,2)+'.'+string(s+2); - } else { - // Remove zeros starting from right end - char* ptr = s + slen - 1; - while (*ptr == '0' && ptr > s) ptr--; - - if (ptr == s) - out = std::string(s, 1); - else - out = std::string(s, 1) + '.' + std::string(s + 1, ptr - (s + 1) + 1); - - // out = string(s,1)+'.'+string(s+1); - } - - // Make final string - if (--exp) { - if (exp > 0) - out += "e+" + mpfr::toString(exp, std::dec); - else - out += "e" + mpfr::toString(exp, std::dec); - } - } - - mpfr_free_str(s); - return out; - } else { - return "conversion error!"; - } -#endif -} - -////////////////////////////////////////////////////////////////////////// -// I/O -inline std::ostream& mpreal::output(std::ostream& os) const { - std::ostringstream format; - const std::ios::fmtflags flags = os.flags(); - - format << ((flags & std::ios::showpos) ? "%+" : "%"); - if (os.precision() >= 0) - format << '.' << os.precision() << "R*" - << ((flags & std::ios::floatfield) == std::ios::fixed - ? 'f' - : (flags & std::ios::floatfield) == std::ios::scientific ? 'e' : 'g'); - else - format << "R*e"; - - char* s = NULL; - if (!(mpfr_asprintf(&s, format.str().c_str(), mpfr::mpreal::get_default_rnd(), mpfr_srcptr()) < 0)) { - os << std::string(s); - mpfr_free_str(s); - } - return os; -} - -inline std::ostream& operator<<(std::ostream& os, const mpreal& v) { return v.output(os); } - -inline std::istream& operator>>(std::istream& is, mpreal& v) { - // TODO: use cout::hexfloat and other flags to setup base - std::string tmp; - is >> tmp; - mpfr_set_str(v.mpfr_ptr(), tmp.c_str(), 10, mpreal::get_default_rnd()); - return is; -} - -////////////////////////////////////////////////////////////////////////// -// Bits - decimal digits relation -// bits = ceil(digits*log[2](10)) -// digits = floor(bits*log[10](2)) - -inline mp_prec_t digits2bits(int d) { - const double LOG2_10 = 3.3219280948873624; - - return mp_prec_t(std::ceil(d * LOG2_10)); -} - -inline int bits2digits(mp_prec_t b) { - const double LOG10_2 = 0.30102999566398119; - - return int(std::floor(b * LOG10_2)); -} - -////////////////////////////////////////////////////////////////////////// -// Set/Get number properties -inline int sgn(const mpreal& op) { return mpfr_sgn(op.mpfr_srcptr()); } - -inline mpreal& mpreal::setSign(int sign, mp_rnd_t RoundingMode) { - mpfr_setsign(mpfr_ptr(), mpfr_srcptr(), (sign < 0 ? 1 : 0), RoundingMode); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline int mpreal::getPrecision() const { return int(mpfr_get_prec(mpfr_srcptr())); } - -inline mpreal& mpreal::setPrecision(int Precision, mp_rnd_t RoundingMode) { - mpfr_prec_round(mpfr_ptr(), Precision, RoundingMode); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::setInf(int sign) { - mpfr_set_inf(mpfr_ptr(), sign); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::setNan() { - mpfr_set_nan(mpfr_ptr()); - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mpreal& mpreal::setZero(int sign) { -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3, 0, 0)) - mpfr_set_zero(mpfr_ptr(), sign); -#else - mpfr_set_si(mpfr_ptr(), 0, (mpfr_get_default_rounding_mode)()); - setSign(sign); -#endif - - MPREAL_MSVC_DEBUGVIEW_CODE; - return *this; -} - -inline mp_prec_t mpreal::get_prec() const { return mpfr_get_prec(mpfr_srcptr()); } - -inline void mpreal::set_prec(mp_prec_t prec, mp_rnd_t rnd_mode) { - mpfr_prec_round(mpfr_ptr(), prec, rnd_mode); - MPREAL_MSVC_DEBUGVIEW_CODE; -} - -inline mp_exp_t mpreal::get_exp() { return mpfr_get_exp(mpfr_srcptr()); } - -inline int mpreal::set_exp(mp_exp_t e) { - int x = mpfr_set_exp(mpfr_ptr(), e); - MPREAL_MSVC_DEBUGVIEW_CODE; - return x; -} - -inline const mpreal frexp(const mpreal& v, mp_exp_t* exp) { - mpreal x(v); - *exp = x.get_exp(); - x.set_exp(0); - return x; -} - -inline const mpreal ldexp(const mpreal& v, mp_exp_t exp) { - mpreal x(v); - - // rounding is not important since we are just increasing the exponent (= exact operation) - mpfr_mul_2si(x.mpfr_ptr(), x.mpfr_srcptr(), exp, mpreal::get_default_rnd()); - return x; -} - -inline const mpreal scalbn(const mpreal& v, mp_exp_t exp) { return ldexp(v, exp); } - -inline mpreal machine_epsilon(mp_prec_t prec) { - /* the smallest eps such that 1 + eps != 1 */ - return machine_epsilon(mpreal(1, prec)); -} - -inline mpreal machine_epsilon(const mpreal& x) { - /* the smallest eps such that x + eps != x */ - if (x < 0) { - return nextabove(-x) + x; - } else { - return nextabove(x) - x; - } -} - -// minval is 'safe' meaning 1 / minval does not overflow -inline mpreal minval(mp_prec_t prec) { - /* min = 1/2 * 2^emin = 2^(emin - 1) */ - return mpreal(1, prec) << mpreal::get_emin() - 1; -} - -// maxval is 'safe' meaning 1 / maxval does not underflow -inline mpreal maxval(mp_prec_t prec) { - /* max = (1 - eps) * 2^emax, eps is machine epsilon */ - return (mpreal(1, prec) - machine_epsilon(prec)) << mpreal::get_emax(); -} - -inline bool isEqualUlps(const mpreal& a, const mpreal& b, int maxUlps) { - return abs(a - b) <= machine_epsilon((max)(abs(a), abs(b))) * maxUlps; -} - -inline bool isEqualFuzzy(const mpreal& a, const mpreal& b, const mpreal& eps) { return abs(a - b) <= eps; } - -inline bool isEqualFuzzy(const mpreal& a, const mpreal& b) { - return isEqualFuzzy(a, b, machine_epsilon((max)(1, (min)(abs(a), abs(b))))); -} - -////////////////////////////////////////////////////////////////////////// -// C++11 sign functions. -inline mpreal copysign(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal rop(0, mpfr_get_prec(x.mpfr_ptr())); - mpfr_setsign(rop.mpfr_ptr(), x.mpfr_srcptr(), mpfr_signbit(y.mpfr_srcptr()), rnd_mode); - return rop; -} - -inline bool signbit(const mpreal& x) { return mpfr_signbit(x.mpfr_srcptr()); } - -inline const mpreal modf(const mpreal& v, mpreal& n) { - mpreal f(v); - - // rounding is not important since we are using the same number - mpfr_frac(f.mpfr_ptr(), f.mpfr_srcptr(), mpreal::get_default_rnd()); - mpfr_trunc(n.mpfr_ptr(), v.mpfr_srcptr()); - return f; -} - -inline int mpreal::check_range(int t, mp_rnd_t rnd_mode) { return mpfr_check_range(mpfr_ptr(), t, rnd_mode); } - -inline int mpreal::subnormalize(int t, mp_rnd_t rnd_mode) { - int r = mpfr_subnormalize(mpfr_ptr(), t, rnd_mode); - MPREAL_MSVC_DEBUGVIEW_CODE; - return r; -} - -inline mp_exp_t mpreal::get_emin(void) { return mpfr_get_emin(); } - -inline int mpreal::set_emin(mp_exp_t exp) { return mpfr_set_emin(exp); } - -inline mp_exp_t mpreal::get_emax(void) { return mpfr_get_emax(); } - -inline int mpreal::set_emax(mp_exp_t exp) { return mpfr_set_emax(exp); } - -inline mp_exp_t mpreal::get_emin_min(void) { return mpfr_get_emin_min(); } - -inline mp_exp_t mpreal::get_emin_max(void) { return mpfr_get_emin_max(); } - -inline mp_exp_t mpreal::get_emax_min(void) { return mpfr_get_emax_min(); } - -inline mp_exp_t mpreal::get_emax_max(void) { return mpfr_get_emax_max(); } - -////////////////////////////////////////////////////////////////////////// -// Mathematical Functions -////////////////////////////////////////////////////////////////////////// -#define MPREAL_UNARY_MATH_FUNCTION_BODY(f) \ - mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); \ - mpfr_##f(y.mpfr_ptr(), x.mpfr_srcptr(), r); \ - return y; - -inline const mpreal sqr(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(sqr); -} - -inline const mpreal sqrt(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(sqrt); -} - -inline const mpreal sqrt(const unsigned long int x, mp_rnd_t r) { - mpreal y; - mpfr_sqrt_ui(y.mpfr_ptr(), x, r); - return y; -} - -inline const mpreal sqrt(const unsigned int v, mp_rnd_t rnd_mode) { - return sqrt(static_cast(v), rnd_mode); -} - -inline const mpreal sqrt(const long int v, mp_rnd_t rnd_mode) { - if (v >= 0) - return sqrt(static_cast(v), rnd_mode); - else - return mpreal().setNan(); // NaN -} - -inline const mpreal sqrt(const int v, mp_rnd_t rnd_mode) { - if (v >= 0) - return sqrt(static_cast(v), rnd_mode); - else - return mpreal().setNan(); // NaN -} - -inline const mpreal root(const mpreal& x, unsigned long int k, mp_rnd_t r = mpreal::get_default_rnd()) { - mpreal y(0, mpfr_get_prec(x.mpfr_srcptr())); - mpfr_rootn_ui(y.mpfr_ptr(), x.mpfr_srcptr(), k, r); - return y; -} - -inline const mpreal dim(const mpreal& a, const mpreal& b, mp_rnd_t r = mpreal::get_default_rnd()) { - mpreal y(0, mpfr_get_prec(a.mpfr_srcptr())); - mpfr_dim(y.mpfr_ptr(), a.mpfr_srcptr(), b.mpfr_srcptr(), r); - return y; -} - -inline int cmpabs(const mpreal& a, const mpreal& b) { return mpfr_cmpabs(a.mpfr_ptr(), b.mpfr_srcptr()); } - -inline int sin_cos(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - return mpfr_sin_cos(s.mpfr_ptr(), c.mpfr_ptr(), v.mpfr_srcptr(), rnd_mode); -} - -inline const mpreal sqrt(const long double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v), rnd_mode); } -inline const mpreal sqrt(const double v, mp_rnd_t rnd_mode) { return sqrt(mpreal(v), rnd_mode); } - -inline const mpreal cbrt(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(cbrt); -} -inline const mpreal fabs(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(abs); -} -inline const mpreal abs(const mpreal& x, mp_rnd_t r) { MPREAL_UNARY_MATH_FUNCTION_BODY(abs); } -inline const mpreal log(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(log); -} -inline const mpreal log2(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(log2); -} -inline const mpreal log10(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(log10); -} -inline const mpreal exp(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(exp); -} -inline const mpreal exp2(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(exp2); -} -inline const mpreal exp10(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(exp10); -} -inline const mpreal cos(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(cos); -} -inline const mpreal sin(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(sin); -} -inline const mpreal tan(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(tan); -} -inline const mpreal sec(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(sec); -} -inline const mpreal csc(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(csc); -} -inline const mpreal cot(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(cot); -} -inline const mpreal acos(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(acos); -} -inline const mpreal asin(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(asin); -} -inline const mpreal atan(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(atan); -} - -inline const mpreal logb(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { return log2(abs(x), r); } - -inline const mpreal acot(const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return atan(1 / v, r); } -inline const mpreal asec(const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return acos(1 / v, r); } -inline const mpreal acsc(const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return asin(1 / v, r); } -inline const mpreal acoth(const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return atanh(1 / v, r); } -inline const mpreal asech(const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return acosh(1 / v, r); } -inline const mpreal acsch(const mpreal& v, mp_rnd_t r = mpreal::get_default_rnd()) { return asinh(1 / v, r); } - -inline const mpreal cosh(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(cosh); -} -inline const mpreal sinh(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(sinh); -} -inline const mpreal tanh(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(tanh); -} -inline const mpreal sech(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(sech); -} -inline const mpreal csch(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(csch); -} -inline const mpreal coth(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(coth); -} -inline const mpreal acosh(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(acosh); -} -inline const mpreal asinh(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(asinh); -} -inline const mpreal atanh(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(atanh); -} - -inline const mpreal log1p(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(log1p); -} -inline const mpreal expm1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(expm1); -} -inline const mpreal eint(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(eint); -} -inline const mpreal gamma(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(gamma); -} -inline const mpreal tgamma(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(gamma); -} -inline const mpreal lngamma(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(lngamma); -} -inline const mpreal zeta(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(zeta); -} -inline const mpreal erf(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(erf); -} -inline const mpreal erfc(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(erfc); -} -inline const mpreal besselj0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(j0); -} -inline const mpreal besselj1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(j1); -} -inline const mpreal bessely0(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(y0); -} -inline const mpreal bessely1(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(y1); -} - -inline const mpreal atan2(const mpreal& y, const mpreal& x, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a(0, (std::max)(y.getPrecision(), x.getPrecision())); - mpfr_atan2(a.mpfr_ptr(), y.mpfr_srcptr(), x.mpfr_srcptr(), rnd_mode); - return a; -} - -inline const mpreal hypot(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a(0, (std::max)(y.getPrecision(), x.getPrecision())); - mpfr_hypot(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode); - return a; -} - -inline const mpreal remainder(const mpreal& x, const mpreal& y, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a(0, (std::max)(y.getPrecision(), x.getPrecision())); - mpfr_remainder(a.mpfr_ptr(), x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode); - return a; -} - -inline const mpreal remquo(long* q, const mpreal& x, const mpreal& y, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a(0, (std::max)(y.getPrecision(), x.getPrecision())); - mpfr_remquo(a.mpfr_ptr(), q, x.mpfr_srcptr(), y.mpfr_srcptr(), rnd_mode); - return a; -} - -inline const mpreal fac_ui(unsigned long int v, mp_prec_t prec = mpreal::get_default_prec(), - mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x(0, prec); - mpfr_fac_ui(x.mpfr_ptr(), v, rnd_mode); - return x; -} - -inline const mpreal lgamma(const mpreal& v, int* signp = 0, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x(v); - int tsignp; - - if (signp) - mpfr_lgamma(x.mpfr_ptr(), signp, v.mpfr_srcptr(), rnd_mode); - else - mpfr_lgamma(x.mpfr_ptr(), &tsignp, v.mpfr_srcptr(), rnd_mode); - - return x; -} - -inline const mpreal besseljn(long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - mpreal y(0, x.getPrecision()); - mpfr_jn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r); - return y; -} - -inline const mpreal besselyn(long n, const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - mpreal y(0, x.getPrecision()); - mpfr_yn(y.mpfr_ptr(), n, x.mpfr_srcptr(), r); - return y; -} - -inline const mpreal fma(const mpreal& v1, const mpreal& v2, const mpreal& v3, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a; - mp_prec_t p1, p2, p3; - - p1 = v1.get_prec(); - p2 = v2.get_prec(); - p3 = v3.get_prec(); - - a.set_prec(p3 > p2 ? (p3 > p1 ? p3 : p1) : (p2 > p1 ? p2 : p1)); - - mpfr_fma(a.mp, v1.mp, v2.mp, v3.mp, rnd_mode); - return a; -} - -inline const mpreal fms(const mpreal& v1, const mpreal& v2, const mpreal& v3, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a; - mp_prec_t p1, p2, p3; - - p1 = v1.get_prec(); - p2 = v2.get_prec(); - p3 = v3.get_prec(); - - a.set_prec(p3 > p2 ? (p3 > p1 ? p3 : p1) : (p2 > p1 ? p2 : p1)); - - mpfr_fms(a.mp, v1.mp, v2.mp, v3.mp, rnd_mode); - return a; -} - -inline const mpreal agm(const mpreal& v1, const mpreal& v2, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a; - mp_prec_t p1, p2; - - p1 = v1.get_prec(); - p2 = v2.get_prec(); - - a.set_prec(p1 > p2 ? p1 : p2); - - mpfr_agm(a.mp, v1.mp, v2.mp, rnd_mode); - - return a; -} - -inline const mpreal sum(const mpreal tab[], const unsigned long int n, int& status, - mp_rnd_t mode = mpreal::get_default_rnd()) { - mpfr_srcptr* p = new mpfr_srcptr[n]; - - for (unsigned long int i = 0; i < n; i++) p[i] = tab[i].mpfr_srcptr(); - - mpreal x; - status = mpfr_sum(x.mpfr_ptr(), (mpfr_ptr*)p, n, mode); - - delete[] p; - return x; -} - -////////////////////////////////////////////////////////////////////////// -// MPFR 2.4.0 Specifics -#if (MPFR_VERSION >= MPFR_VERSION_NUM(2, 4, 0)) - -inline int sinh_cosh(mpreal& s, mpreal& c, const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - return mpfr_sinh_cosh(s.mp, c.mp, v.mp, rnd_mode); -} - -inline const mpreal li2(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(li2); -} - -inline const mpreal rem(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - /* R = rem(X,Y) if Y != 0, returns X - n * Y where n = trunc(X/Y). */ - return fmod(x, y, rnd_mode); -} - -inline const mpreal mod(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - (void)rnd_mode; - - /* - - m = mod(x,y) if y != 0, returns x - n*y where n = floor(x/y) - - The following are true by convention: - - mod(x,0) is x - - mod(x,x) is 0 - - mod(x,y) for x != y and y != 0 has the same sign as y. - - */ - - if (iszero(y)) return x; - if (x == y) return 0; - - mpreal m = x - floor(x / y) * y; - - m.setSign(sgn(y)); // make sure result has the same sign as Y - - return m; -} - -inline const mpreal fmod(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a; - mp_prec_t yp, xp; - - yp = y.get_prec(); - xp = x.get_prec(); - - a.set_prec(yp > xp ? yp : xp); - - mpfr_fmod(a.mp, x.mp, y.mp, rnd_mode); - - return a; -} - -inline const mpreal rec_sqrt(const mpreal& v, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x(v); - mpfr_rec_sqrt(x.mp, v.mp, rnd_mode); - return x; -} -#endif // MPFR 2.4.0 Specifics - -////////////////////////////////////////////////////////////////////////// -// MPFR 3.0.0 Specifics -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3, 0, 0)) -inline const mpreal digamma(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(digamma); -} -inline const mpreal ai(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(ai); -} -#endif // MPFR 3.0.0 Specifics - -////////////////////////////////////////////////////////////////////////// -// Constants -inline const mpreal const_log2(mp_prec_t p = mpreal::get_default_prec(), - mp_rnd_t r = mpreal::get_default_rnd()) { - mpreal x(0, p); - mpfr_const_log2(x.mpfr_ptr(), r); - return x; -} - -inline const mpreal const_pi(mp_prec_t p = mpreal::get_default_prec(), - mp_rnd_t r = mpreal::get_default_rnd()) { - mpreal x(0, p); - mpfr_const_pi(x.mpfr_ptr(), r); - return x; -} - -inline const mpreal const_euler(mp_prec_t p = mpreal::get_default_prec(), - mp_rnd_t r = mpreal::get_default_rnd()) { - mpreal x(0, p); - mpfr_const_euler(x.mpfr_ptr(), r); - return x; -} - -inline const mpreal const_catalan(mp_prec_t p = mpreal::get_default_prec(), - mp_rnd_t r = mpreal::get_default_rnd()) { - mpreal x(0, p); - mpfr_const_catalan(x.mpfr_ptr(), r); - return x; -} - -inline const mpreal const_infinity(int sign = 1, mp_prec_t p = mpreal::get_default_prec()) { - mpreal x(0, p); - mpfr_set_inf(x.mpfr_ptr(), sign); - return x; -} - -////////////////////////////////////////////////////////////////////////// -// Integer Related Functions -inline const mpreal ceil(const mpreal& v) { - mpreal x(v); - mpfr_ceil(x.mp, v.mp); - return x; -} - -inline const mpreal floor(const mpreal& v) { - mpreal x(v); - mpfr_floor(x.mp, v.mp); - return x; -} - -inline const mpreal round(const mpreal& v) { - mpreal x(v); - mpfr_round(x.mp, v.mp); - return x; -} - -inline const mpreal trunc(const mpreal& v) { - mpreal x(v); - mpfr_trunc(x.mp, v.mp); - return x; -} - -inline const mpreal rint(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(rint); -} -inline const mpreal rint_ceil(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(rint_ceil); -} -inline const mpreal rint_floor(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(rint_floor); -} -inline const mpreal rint_round(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(rint_round); -} -inline const mpreal rint_trunc(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(rint_trunc); -} -inline const mpreal frac(const mpreal& x, mp_rnd_t r = mpreal::get_default_rnd()) { - MPREAL_UNARY_MATH_FUNCTION_BODY(frac); -} - -////////////////////////////////////////////////////////////////////////// -// Miscellaneous Functions -inline void swap(mpreal& a, mpreal& b) { mpfr_swap(a.mp, b.mp); } -inline const mpreal(max)(const mpreal& x, const mpreal& y) { return (x > y ? x : y); } -inline const mpreal(min)(const mpreal& x, const mpreal& y) { return (x < y ? x : y); } - -inline const mpreal fmax(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a; - mpfr_max(a.mp, x.mp, y.mp, rnd_mode); - return a; -} - -inline const mpreal fmin(const mpreal& x, const mpreal& y, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal a; - mpfr_min(a.mp, x.mp, y.mp, rnd_mode); - return a; -} - -inline const mpreal nexttoward(const mpreal& x, const mpreal& y) { - mpreal a(x); - mpfr_nexttoward(a.mp, y.mp); - return a; -} - -inline const mpreal nextabove(const mpreal& x) { - mpreal a(x); - mpfr_nextabove(a.mp); - return a; -} - -inline const mpreal nextbelow(const mpreal& x) { - mpreal a(x); - mpfr_nextbelow(a.mp); - return a; -} - -inline const mpreal urandomb(gmp_randstate_t& state) { - mpreal x; - mpfr_urandomb(x.mpfr_ptr(), state); - return x; -} - -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3, 0, 0)) -inline const mpreal urandom(gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x; - mpfr_urandom(x.mpfr_ptr(), state, rnd_mode); - return x; -} -#endif - -#if (MPFR_VERSION <= MPFR_VERSION_NUM(2, 4, 2)) -inline const mpreal random2(mp_size_t size, mp_exp_t exp) { - mpreal x; - mpfr_random2(x.mpfr_ptr(), size, exp); - return x; -} -#endif - -// Uniformly distributed random number generation -// a = random(seed); <- initialization & first random number generation -// a = random(); <- next random numbers generation -// seed != 0 -inline const mpreal random(unsigned int seed = 0) { -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3, 0, 0)) - static gmp_randstate_t state; - static bool initialize = true; - - if (initialize) { - gmp_randinit_default(state); - gmp_randseed_ui(state, 0); - initialize = false; - } - - if (seed != 0) gmp_randseed_ui(state, seed); - - return mpfr::urandom(state); -#else - if (seed != 0) std::srand(seed); - return mpfr::mpreal(std::rand() / (double)RAND_MAX); -#endif -} - -#if (MPFR_VERSION >= MPFR_VERSION_NUM(3, 1, 0)) - -inline const mpreal grandom(gmp_randstate_t& state, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x; - mpfr_nrandom(x.mpfr_ptr(), state, rnd_mode); - return x; -} - -inline const mpreal grandom(unsigned int seed = 0) { - static gmp_randstate_t state; - static bool initialize = true; - - if (initialize) { - gmp_randinit_default(state); - gmp_randseed_ui(state, 0); - initialize = false; - } - - if (seed != 0) gmp_randseed_ui(state, seed); - - return mpfr::grandom(state); -} -#endif - -////////////////////////////////////////////////////////////////////////// -// Set/Get global properties -inline void mpreal::set_default_prec(mp_prec_t prec) { mpfr_set_default_prec(prec); } - -inline void mpreal::set_default_rnd(mp_rnd_t rnd_mode) { mpfr_set_default_rounding_mode(rnd_mode); } - -inline bool mpreal::fits_in_bits(double x, int n) { - int i; - double t; - return IsInf(x) || (std::modf(std::ldexp(std::frexp(x, &i), n), &t) == 0.0); -} - -inline const mpreal pow(const mpreal& a, const mpreal& b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x(a); - mpfr_pow(x.mp, x.mp, b.mp, rnd_mode); - return x; -} - -inline const mpreal pow(const mpreal& a, const mpz_t b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x(a); - mpfr_pow_z(x.mp, x.mp, b, rnd_mode); - return x; -} - -inline const mpreal pow(const mpreal& a, const unsigned long int b, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x(a); - mpfr_pow_ui(x.mp, x.mp, b, rnd_mode); - return x; -} - -inline const mpreal pow(const mpreal& a, const unsigned int b, mp_rnd_t rnd_mode) { - return pow(a, static_cast(b), rnd_mode); -} - -inline const mpreal pow(const mpreal& a, const long int b, mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x(a); - mpfr_pow_si(x.mp, x.mp, b, rnd_mode); - return x; -} - -inline const mpreal pow(const mpreal& a, const int b, mp_rnd_t rnd_mode) { - return pow(a, static_cast(b), rnd_mode); -} - -inline const mpreal pow(const mpreal& a, const long double b, mp_rnd_t rnd_mode) { - return pow(a, mpreal(b), rnd_mode); -} - -inline const mpreal pow(const mpreal& a, const double b, mp_rnd_t rnd_mode) { - return pow(a, mpreal(b), rnd_mode); -} - -inline const mpreal pow(const unsigned long int a, const mpreal& b, - mp_rnd_t rnd_mode = mpreal::get_default_rnd()) { - mpreal x(a); - mpfr_ui_pow(x.mp, a, b.mp, rnd_mode); - return x; -} - -inline const mpreal pow(const unsigned int a, const mpreal& b, mp_rnd_t rnd_mode) { - return pow(static_cast(a), b, rnd_mode); -} - -inline const mpreal pow(const long int a, const mpreal& b, mp_rnd_t rnd_mode) { - if (a >= 0) - return pow(static_cast(a), b, rnd_mode); - else - return pow(mpreal(a), b, rnd_mode); -} - -inline const mpreal pow(const int a, const mpreal& b, mp_rnd_t rnd_mode) { - if (a >= 0) - return pow(static_cast(a), b, rnd_mode); - else - return pow(mpreal(a), b, rnd_mode); -} - -inline const mpreal pow(const long double a, const mpreal& b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), b, rnd_mode); -} - -inline const mpreal pow(const double a, const mpreal& b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), b, rnd_mode); -} - -// pow unsigned long int -inline const mpreal pow(const unsigned long int a, const unsigned long int b, mp_rnd_t rnd_mode) { - mpreal x(a); - mpfr_ui_pow_ui(x.mp, a, b, rnd_mode); - return x; -} - -inline const mpreal pow(const unsigned long int a, const unsigned int b, mp_rnd_t rnd_mode) { - return pow(a, static_cast(b), rnd_mode); // mpfr_ui_pow_ui -} - -inline const mpreal pow(const unsigned long int a, const long int b, mp_rnd_t rnd_mode) { - if (b > 0) - return pow(a, static_cast(b), rnd_mode); // mpfr_ui_pow_ui - else - return pow(a, mpreal(b), rnd_mode); // mpfr_ui_pow -} - -inline const mpreal pow(const unsigned long int a, const int b, mp_rnd_t rnd_mode) { - if (b > 0) - return pow(a, static_cast(b), rnd_mode); // mpfr_ui_pow_ui - else - return pow(a, mpreal(b), rnd_mode); // mpfr_ui_pow -} - -inline const mpreal pow(const unsigned long int a, const long double b, mp_rnd_t rnd_mode) { - return pow(a, mpreal(b), rnd_mode); // mpfr_ui_pow -} - -inline const mpreal pow(const unsigned long int a, const double b, mp_rnd_t rnd_mode) { - return pow(a, mpreal(b), rnd_mode); // mpfr_ui_pow -} - -// pow unsigned int -inline const mpreal pow(const unsigned int a, const unsigned long int b, mp_rnd_t rnd_mode) { - return pow(static_cast(a), b, rnd_mode); // mpfr_ui_pow_ui -} - -inline const mpreal pow(const unsigned int a, const unsigned int b, mp_rnd_t rnd_mode) { - return pow(static_cast(a), static_cast(b), - rnd_mode); // mpfr_ui_pow_ui -} - -inline const mpreal pow(const unsigned int a, const long int b, mp_rnd_t rnd_mode) { - if (b > 0) - return pow(static_cast(a), static_cast(b), - rnd_mode); // mpfr_ui_pow_ui - else - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow -} - -inline const mpreal pow(const unsigned int a, const int b, mp_rnd_t rnd_mode) { - if (b > 0) - return pow(static_cast(a), static_cast(b), - rnd_mode); // mpfr_ui_pow_ui - else - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow -} - -inline const mpreal pow(const unsigned int a, const long double b, mp_rnd_t rnd_mode) { - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow -} - -inline const mpreal pow(const unsigned int a, const double b, mp_rnd_t rnd_mode) { - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow -} - -// pow long int -inline const mpreal pow(const long int a, const unsigned long int b, mp_rnd_t rnd_mode) { - if (a > 0) - return pow(static_cast(a), b, rnd_mode); // mpfr_ui_pow_ui - else - return pow(mpreal(a), b, rnd_mode); // mpfr_pow_ui -} - -inline const mpreal pow(const long int a, const unsigned int b, mp_rnd_t rnd_mode) { - if (a > 0) - return pow(static_cast(a), static_cast(b), - rnd_mode); // mpfr_ui_pow_ui - else - return pow(mpreal(a), static_cast(b), rnd_mode); // mpfr_pow_ui -} - -inline const mpreal pow(const long int a, const long int b, mp_rnd_t rnd_mode) { - if (a > 0) { - if (b > 0) - return pow(static_cast(a), static_cast(b), - rnd_mode); // mpfr_ui_pow_ui - else - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow - } else { - return pow(mpreal(a), b, rnd_mode); // mpfr_pow_si - } -} - -inline const mpreal pow(const long int a, const int b, mp_rnd_t rnd_mode) { - if (a > 0) { - if (b > 0) - return pow(static_cast(a), static_cast(b), - rnd_mode); // mpfr_ui_pow_ui - else - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow - } else { - return pow(mpreal(a), static_cast(b), rnd_mode); // mpfr_pow_si - } -} - -inline const mpreal pow(const long int a, const long double b, mp_rnd_t rnd_mode) { - if (a >= 0) - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow - else - return pow(mpreal(a), mpreal(b), rnd_mode); // mpfr_pow -} - -inline const mpreal pow(const long int a, const double b, mp_rnd_t rnd_mode) { - if (a >= 0) - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow - else - return pow(mpreal(a), mpreal(b), rnd_mode); // mpfr_pow -} - -// pow int -inline const mpreal pow(const int a, const unsigned long int b, mp_rnd_t rnd_mode) { - if (a > 0) - return pow(static_cast(a), b, rnd_mode); // mpfr_ui_pow_ui - else - return pow(mpreal(a), b, rnd_mode); // mpfr_pow_ui -} - -inline const mpreal pow(const int a, const unsigned int b, mp_rnd_t rnd_mode) { - if (a > 0) - return pow(static_cast(a), static_cast(b), - rnd_mode); // mpfr_ui_pow_ui - else - return pow(mpreal(a), static_cast(b), rnd_mode); // mpfr_pow_ui -} - -inline const mpreal pow(const int a, const long int b, mp_rnd_t rnd_mode) { - if (a > 0) { - if (b > 0) - return pow(static_cast(a), static_cast(b), - rnd_mode); // mpfr_ui_pow_ui - else - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow - } else { - return pow(mpreal(a), b, rnd_mode); // mpfr_pow_si - } -} - -inline const mpreal pow(const int a, const int b, mp_rnd_t rnd_mode) { - if (a > 0) { - if (b > 0) - return pow(static_cast(a), static_cast(b), - rnd_mode); // mpfr_ui_pow_ui - else - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow - } else { - return pow(mpreal(a), static_cast(b), rnd_mode); // mpfr_pow_si - } -} - -inline const mpreal pow(const int a, const long double b, mp_rnd_t rnd_mode) { - if (a >= 0) - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow - else - return pow(mpreal(a), mpreal(b), rnd_mode); // mpfr_pow -} - -inline const mpreal pow(const int a, const double b, mp_rnd_t rnd_mode) { - if (a >= 0) - return pow(static_cast(a), mpreal(b), rnd_mode); // mpfr_ui_pow - else - return pow(mpreal(a), mpreal(b), rnd_mode); // mpfr_pow -} - -// pow long double -inline const mpreal pow(const long double a, const long double b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), mpreal(b), rnd_mode); -} - -inline const mpreal pow(const long double a, const unsigned long int b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), b, rnd_mode); // mpfr_pow_ui -} - -inline const mpreal pow(const long double a, const unsigned int b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), static_cast(b), rnd_mode); // mpfr_pow_ui -} - -inline const mpreal pow(const long double a, const long int b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), b, rnd_mode); // mpfr_pow_si -} - -inline const mpreal pow(const long double a, const int b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), static_cast(b), rnd_mode); // mpfr_pow_si -} - -inline const mpreal pow(const double a, const double b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), mpreal(b), rnd_mode); -} - -inline const mpreal pow(const double a, const unsigned long int b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), b, rnd_mode); // mpfr_pow_ui -} - -inline const mpreal pow(const double a, const unsigned int b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), static_cast(b), rnd_mode); // mpfr_pow_ui -} - -inline const mpreal pow(const double a, const long int b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), b, rnd_mode); // mpfr_pow_si -} - -inline const mpreal pow(const double a, const int b, mp_rnd_t rnd_mode) { - return pow(mpreal(a), static_cast(b), rnd_mode); // mpfr_pow_si -} -} // namespace mpfr - -// Explicit specialization of std::swap for mpreal numbers -// Thus standard algorithms will use efficient version of swap (due to Koenig lookup) -// Non-throwing swap C++ idiom: http://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Non-throwing_swap -namespace std { -// we are allowed to extend namespace std with specializations only -template <> -inline void swap(mpfr::mpreal& x, mpfr::mpreal& y) { - return mpfr::swap(x, y); -} - -template <> -class numeric_limits { - public: - static const bool is_specialized = true; - static const bool is_signed = true; - static const bool is_integer = false; - static const bool is_exact = false; - static const int radix = 2; - - static const bool has_infinity = true; - static const bool has_quiet_NaN = true; - static const bool has_signaling_NaN = true; - - static const bool is_iec559 = true; // = IEEE 754 - static const bool is_bounded = true; - static const bool is_modulo = false; - static const bool traps = true; - static const bool tinyness_before = true; - - static const float_denorm_style has_denorm = denorm_absent; - - inline static mpfr::mpreal(min)(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { - return mpfr::minval(precision); - } - inline static mpfr::mpreal(max)(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { - return mpfr::maxval(precision); - } - inline static mpfr::mpreal lowest(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { - return -mpfr::maxval(precision); - } - - // Returns smallest eps such that 1 + eps != 1 (classic machine epsilon) - inline static mpfr::mpreal epsilon(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { - return mpfr::machine_epsilon(precision); - } - - // Returns smallest eps such that x + eps != x (relative machine epsilon) - inline static mpfr::mpreal epsilon(const mpfr::mpreal& x) { return mpfr::machine_epsilon(x); } - - inline static mpfr::mpreal round_error(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { - mp_rnd_t r = mpfr::mpreal::get_default_rnd(); - - if (r == GMP_RNDN) - return mpfr::mpreal(0.5, precision); - else - return mpfr::mpreal(1.0, precision); - } - - inline static const mpfr::mpreal infinity() { return mpfr::const_infinity(); } - inline static const mpfr::mpreal quiet_NaN() { return mpfr::mpreal().setNan(); } - inline static const mpfr::mpreal signaling_NaN() { return mpfr::mpreal().setNan(); } - inline static const mpfr::mpreal denorm_min() { return (min)(); } - - // Please note, exponent range is not fixed in MPFR - static const int min_exponent = MPFR_EMIN_DEFAULT; - static const int max_exponent = MPFR_EMAX_DEFAULT; - MPREAL_PERMISSIVE_EXPR static const int min_exponent10 = (int)(MPFR_EMIN_DEFAULT * 0.3010299956639811); - MPREAL_PERMISSIVE_EXPR static const int max_exponent10 = (int)(MPFR_EMAX_DEFAULT * 0.3010299956639811); - -#ifdef MPREAL_HAVE_DYNAMIC_STD_NUMERIC_LIMITS - - // Following members should be constant according to standard, but they can be variable in MPFR - // So we define them as functions here. - // - // This is preferable way for std::numeric_limits specialization. - // But it is incompatible with standard std::numeric_limits and might not work with other libraries, e.g. - // boost. - // See below for compatible implementation. - inline static float_round_style round_style() { - mp_rnd_t r = mpfr::mpreal::get_default_rnd(); - - switch (r) { - case GMP_RNDN: - return round_to_nearest; - case GMP_RNDZ: - return round_toward_zero; - case GMP_RNDU: - return round_toward_infinity; - case GMP_RNDD: - return round_toward_neg_infinity; - default: - return round_indeterminate; - } - } - - inline static int digits() { return int(mpfr::mpreal::get_default_prec()); } - inline static int digits(const mpfr::mpreal& x) { return x.getPrecision(); } - - inline static int digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { - return mpfr::bits2digits(precision); - } - - inline static int digits10(const mpfr::mpreal& x) { return mpfr::bits2digits(x.getPrecision()); } - - inline static int max_digits10(mp_prec_t precision = mpfr::mpreal::get_default_prec()) { - return digits10(precision); - } -#else - // Digits and round_style are NOT constants when it comes to mpreal. - // If possible, please use functions digits() and round_style() defined above. - // - // These (default) values are preserved for compatibility with existing libraries, e.g. boost. - // Change them accordingly to your application. - // - // For example, if you use 256 bits of precision uniformly in your program, then: - // digits = 256 - // digits10 = 77 - // max_digits10 = 78 - // - // Approximate formula for decimal digits is: digits10 = floor(log10(2) * digits). See bits2digits() for - // more details. - - static const std::float_round_style round_style = round_to_nearest; - static const int digits = 53; - static const int digits10 = 15; - static const int max_digits10 = 16; -#endif -}; -} // namespace std - -#endif /* __MPREAL_H__ */ diff --git a/pipeline/axionGenerator.rml b/pipeline/axionGenerator.rml deleted file mode 100644 index 037c81cf..00000000 --- a/pipeline/axionGenerator.rml +++ /dev/null @@ -1,113 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/metadata/magneticField/magneticField.py b/pipeline/metadata/magneticField/magneticField.py index 1e299d7e..199d4436 100755 --- a/pipeline/metadata/magneticField/magneticField.py +++ b/pipeline/metadata/magneticField/magneticField.py @@ -33,6 +33,16 @@ b3 = int( 1000 * babyField.GetTransversalComponent( p3, d1)) b4 = int( 1000 * babyField.GetTransversalComponent( p4, d1)) +startFrom = ROOT.TVector3(-300,300,-2500) +goTo = ROOT.TVector3(300,-300,2500) +b5 = int( 1000. * babyField.GetTransversalFieldAverage( startFrom, goTo, 20, 500) ) + +print ( "\nEvaluating transversal field average Bykovskiy_201906.dat centered at (0,0,0)" ) +if( b5 != 1548 ): + print ("\nEvaluation of field failed! Exit code : 302") + exit(302) +print ( "[\033[92m OK \x1b[0m]" ) + print ( "\nEvaluating field volume Bykovskiy_201906.dat centered at (0,0,0)" ) if( b1 != 2007 or b2 != 1998 or b3 != 1998 or b4 != 1119 ): print ("\nEvaluation of field failed! Exit code : 102") diff --git a/pipeline/metadata/optics/basic.py b/pipeline/metadata/optics/basic.py index ce0006f3..ff17bc7d 100755 --- a/pipeline/metadata/optics/basic.py +++ b/pipeline/metadata/optics/basic.py @@ -1,5 +1,7 @@ #!/usr/bin/python3 +### This script is now obsolete. It should be updated to use a specific optics inherited class. + import ROOT, math outfname = "opticsBasic.png" diff --git a/pipeline/metadata/optics/optics.py b/pipeline/metadata/optics/mcpl.py similarity index 100% rename from pipeline/metadata/optics/optics.py rename to pipeline/metadata/optics/mcpl.py diff --git a/pipeline/ray-tracing/ValidateXMM.C b/pipeline/ray-tracing/ValidateXMM.C new file mode 100644 index 00000000..2cca73f4 --- /dev/null +++ b/pipeline/ray-tracing/ValidateXMM.C @@ -0,0 +1,33 @@ +#include +#include + +Int_t ValidateXMM(std::string fname) { + std::cout << "Filename : " << fname << std::endl; + TRestRun* run = new TRestRun(fname); + + if (run->GetEntries() != 1000) { + std::cout << "Error. Number of entries is not 1000!" << std::endl; + + return 1; + } + + run->GetAnalysisTree()->Draw("optics_efficiency", "optics_efficiency"); + + 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 << "Efficiency : " << integral / run->GetEntries() << std::endl; + if (integral < 150 || integral > 250) { + std::cout << "Optics efficiency is not within the expected range!" << std::endl; + return 3; + } + + delete run; + + return 0; +} diff --git a/src/TRestAxionAnalysisProcess.cxx b/src/TRestAxionAnalysisProcess.cxx index d1d691e6..4a53e24d 100644 --- a/src/TRestAxionAnalysisProcess.cxx +++ b/src/TRestAxionAnalysisProcess.cxx @@ -114,18 +114,14 @@ TRestEvent* TRestAxionAnalysisProcess::ProcessEvent(TRestEvent* evInput) { SetObservableValue("energy", fAxionEvent->GetEnergy()); - SetObservableValue("posX", fAxionEvent->GetPosition()->X()); - SetObservableValue("posY", fAxionEvent->GetPosition()->Y()); - SetObservableValue("posZ", fAxionEvent->GetPosition()->Z()); + SetObservableValue("posX", fAxionEvent->GetPosition().X()); + SetObservableValue("posY", fAxionEvent->GetPosition().Y()); + SetObservableValue("posZ", fAxionEvent->GetPosition().Z()); - SetObservableValue("probability", fAxionEvent->GetGammaProbability()); + // SetObservableValue("B2", fAxionEvent->GetBSquared()); + // SetObservableValue("Lcoh", fAxionEvent->GetLConversion()); if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fAxionEvent->PrintEvent(); return fAxionEvent; } - -/////////////////////////////////////////////// -/// \brief Function reading input parameters from the RML TRestAxionAnalysisProcess metadata section -/// -void TRestAxionAnalysisProcess::InitFromConfigFile() {} diff --git a/src/TRestAxionOpticsResponseProcess.cxx b/src/TRestAxionDeviationProcess.cxx similarity index 69% rename from src/TRestAxionOpticsResponseProcess.cxx rename to src/TRestAxionDeviationProcess.cxx index 3836dc5a..fdb93cb8 100644 --- a/src/TRestAxionOpticsResponseProcess.cxx +++ b/src/TRestAxionDeviationProcess.cxx @@ -21,7 +21,13 @@ *************************************************************************/ ////////////////////////////////////////////////////////////////////////// -/// TRestAxionOpticsResponseProcess TOBE documented +/// TRestAxionDeviationProcess simply produces a deviation given by an +/// angle, fDevAngle, that defines the cone directrix delimiting the +/// deviation along the main axion direction. +/// +/// For the moment we simply produce a uniform distribution for the angle. +/// This will likely produce a non-homogeneous angular distribution, but +/// still good for small angles. /// ///-------------------------------------------------------------------------- /// @@ -29,23 +35,23 @@ /// /// History of developments: /// -/// 2019-March: First implementation of a dummy optics response +/// 2022-July: Basic implementation of a photon transport process /// Javier Galan /// -/// \class TRestAxionOpticsResponseProcess +/// \class TRestAxionDeviationProcess /// \author /// ///
/// -#include "TRestAxionOpticsResponseProcess.h" +#include "TRestAxionDeviationProcess.h" using namespace std; -ClassImp(TRestAxionOpticsResponseProcess); +ClassImp(TRestAxionDeviationProcess); /////////////////////////////////////////////// /// \brief Default constructor /// -TRestAxionOpticsResponseProcess::TRestAxionOpticsResponseProcess() { Initialize(); } +TRestAxionDeviationProcess::TRestAxionDeviationProcess() { Initialize(); } /////////////////////////////////////////////// /// \brief Constructor loading data from a config file @@ -58,7 +64,7 @@ TRestAxionOpticsResponseProcess::TRestAxionOpticsResponseProcess() { Initialize( /// /// \param cfgFileName A const char* giving the path to an RML file. /// -TRestAxionOpticsResponseProcess::TRestAxionOpticsResponseProcess(char* cfgFileName) { +TRestAxionDeviationProcess::TRestAxionDeviationProcess(char* cfgFileName) { Initialize(); LoadConfig(cfgFileName); @@ -67,12 +73,12 @@ TRestAxionOpticsResponseProcess::TRestAxionOpticsResponseProcess(char* cfgFileNa /////////////////////////////////////////////// /// \brief Default destructor /// -TRestAxionOpticsResponseProcess::~TRestAxionOpticsResponseProcess() { delete fAxionEvent; } +TRestAxionDeviationProcess::~TRestAxionDeviationProcess() { delete fAxionEvent; } /////////////////////////////////////////////// /// \brief Function to load the default config in absence of RML input /// -void TRestAxionOpticsResponseProcess::LoadDefaultConfig() { +void TRestAxionDeviationProcess::LoadDefaultConfig() { SetName(this->ClassName()); SetTitle("Default config"); } @@ -87,26 +93,56 @@ void TRestAxionOpticsResponseProcess::LoadDefaultConfig() { /// \param name The name of the specific metadata. It will be used to find the /// correspondig TRestGeant4AnalysisProcess section inside the RML. /// -void TRestAxionOpticsResponseProcess::LoadConfig(std::string cfgFilename, std::string name) { +void TRestAxionDeviationProcess::LoadConfig(std::string cfgFilename, std::string name) { if (LoadConfigFromFile(cfgFilename, name)) LoadDefaultConfig(); } /////////////////////////////////////////////// /// \brief Function to initialize input/output event members and define the section name /// -void TRestAxionOpticsResponseProcess::Initialize() { +void TRestAxionDeviationProcess::Initialize() { SetSectionName(this->ClassName()); SetLibraryVersion(LIBRARY_VERSION); fAxionEvent = new TRestAxionEvent(); } +/////////////////////////////////////////////// +/// \brief Process initialization. Data members that require initialization just before start processing +/// should be initialized here. +/// +void TRestAxionDeviationProcess::InitProcess() { + RESTDebug << "Entering ... TRestAxionGeneratorProcess::InitProcess" << RESTendl; + + if (!fRandom) { + delete fRandom; + fRandom = nullptr; + } + + fRandom = new TRandom3(fSeed); + if (fSeed == 0) fSeed = fRandom->GetSeed(); +} + /////////////////////////////////////////////// /// \brief The main processing event function /// -TRestEvent* TRestAxionOpticsResponseProcess::ProcessEvent(TRestEvent* evInput) { +TRestEvent* TRestAxionDeviationProcess::ProcessEvent(TRestEvent* evInput) { fAxionEvent = (TRestAxionEvent*)evInput; + TVector3 inPos = fAxionEvent->GetPosition(); + TVector3 inDir = fAxionEvent->GetDirection(); + + Double_t theta = fDevAngle * fRandom->Rndm(); + Double_t phi = TMath::Pi() * fRandom->Rndm(); + + TVector3 ortho = inDir.Orthogonal(); + TVector3 newDir = inDir; + + newDir.Rotate(theta, ortho); + newDir.Rotate(phi, inDir); + + fAxionEvent->SetDirection(newDir); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { fAxionEvent->PrintEvent(); @@ -115,8 +151,3 @@ TRestEvent* TRestAxionOpticsResponseProcess::ProcessEvent(TRestEvent* evInput) { return fAxionEvent; } - -/////////////////////////////////////////////// -/// \brief Function reading input parameters from the RML TRestAxionOpticsResponseProcess metadata section -/// -void TRestAxionOpticsResponseProcess::InitFromConfigFile() {} diff --git a/src/TRestAxionEvent.cxx b/src/TRestAxionEvent.cxx index 64d7a53b..71ef12e8 100644 --- a/src/TRestAxionEvent.cxx +++ b/src/TRestAxionEvent.cxx @@ -82,6 +82,79 @@ TPad* TRestAxionEvent::DrawEvent(TString option) { return fPad; } +/////////////////////////////////////////////// +/// \brief This method will produce a rotation respect to a `center` given by argument. +/// First we rotate the particle and direction along the X-axis by an angle `theta`, then +/// we rotate the particle and direction along the Z-axis by an angle `phi`. +/// +void TRestAxionEvent::RotateXZ(const TVector3& center, Double_t theta, Double_t phi) { + TVector3 ref = fPosition - center; + + ref.RotateX(theta); + ref.RotateZ(phi); + + fPosition = ref + center; + + fDirection.RotateX(theta); + fDirection.RotateZ(phi); +} + +/////////////////////////////////////////////// +/// \brief This method will produce a rotation respect to a `center` given by argument. +/// First we rotate the particle and direction along the Z-axis by an angle `phi`, then +/// we rotate the particle and direction along the X-axis by an angle `theta`. +/// +void TRestAxionEvent::RotateZX(const TVector3& center, Double_t phi, Double_t theta) { + TVector3 ref = fPosition - center; + + ref.RotateZ(phi); + ref.RotateX(theta); + + fPosition = ref + center; + + fDirection.RotateZ(phi); + fDirection.RotateX(theta); +} + +/////////////////////////////////////////////// +/// \brief This method will produce a rotation respect to a `center` given by argument. +/// First we rotate the particle and direction along the X-axis by an angle `pitch`, then +/// we rotate the particle and direction along the Y-axis by an angle `yaw`. +/// +void TRestAxionEvent::RotateXY(const TVector3& center, Double_t pitch, Double_t yaw) { + TVector3 ref = fPosition - center; + + ref.RotateX(pitch); + ref.RotateY(yaw); + + fPosition = ref + center; + + fDirection.RotateX(pitch); + fDirection.RotateY(yaw); +} + +/////////////////////////////////////////////// +/// \brief This method will produce a rotation respect to a `center` given by argument. +/// First we rotate the particle and direction along the Y-axis by an angle `yaw`, then +/// we rotate the particle and direction along the X-axis by an angle `pitch`. +/// +void TRestAxionEvent::RotateYX(const TVector3& center, Double_t yaw, Double_t pitch) { + TVector3 ref = fPosition - center; + + ref.RotateY(yaw); + ref.RotateX(pitch); + + fPosition = ref + center; + + fDirection.RotateY(yaw); + fDirection.RotateX(pitch); +} + +/////////////////////////////////////////////// +/// \brief This method will produce a tranlation of the axion position by an amount `delta`. +/// +void TRestAxionEvent::Translate(const TVector3& delta) { fPosition += delta; } + void TRestAxionEvent::PrintEvent() { TRestEvent::PrintEvent(); @@ -90,6 +163,6 @@ void TRestAxionEvent::PrintEvent() { << endl; cout << "Direction : ( " << fDirection.X() << ", " << fDirection.Y() << ", " << fDirection.Z() << " )" << endl; - cout << "Gamma state probability : " << fGammaProbability << endl; + cout << "B^2 : " << fBSquared << " T^2" << endl; cout << endl; } diff --git a/src/TRestAxionEventProcess.cxx b/src/TRestAxionEventProcess.cxx new file mode 100644 index 00000000..b303a12a --- /dev/null +++ b/src/TRestAxionEventProcess.cxx @@ -0,0 +1,131 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// +/// TOBE written +/// +/// \class TRestAxionEventProcess +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2022-March: First concept. +/// Javier Galan +/// +///
+////////////////////////////////////////////////////////////////////////// + +#include "TRestAxionEventProcess.h" + +using namespace std; + +ClassImp(TRestAxionEventProcess); + +////////////////////////////////////////////////////////////////////////// +/// TRestAxionEventProcess default constructor +/// +TRestAxionEventProcess::TRestAxionEventProcess() { fSingleThreadOnly = false; } + +////////////////////////////////////////////////////////////////////////// +/// TRestAxionEventProcess destructor +/// +TRestAxionEventProcess::~TRestAxionEventProcess() {} + +////////////////////////////////////////////////////////////////////////// +/// \brief Begin of event process, preparation work. Called right before ProcessEvent() +/// +/// This method is called before calling ProcessEvent(). We initialize the process's output +/// event if not null and not same as input event. The event's basic info (ID, timestamp, etc.) +/// will also be set to the same as input event +/// +void TRestAxionEventProcess::BeginOfEventProcess(TRestEvent* inEv) { + TRestEventProcess::BeginOfEventProcess(inEv); + + fAxionEvent = (TRestAxionEvent*)inEv; + RESTDebug << "BoEP: Initial Position. X: " << fAxionEvent->GetPosition().X() + << " Y: " << fAxionEvent->GetPosition().Y() << " Z: " << fAxionEvent->GetPosition().Z() + << RESTendl; + RESTDebug << "BoEP: Initial Direction. X: " << fAxionEvent->GetDirection().X() + << " Y: " << fAxionEvent->GetDirection().Y() << " Z: " << fAxionEvent->GetDirection().Z() + << RESTendl; + fAxionEvent->RotateYX(fCenter, -fYaw, -fPitch); + fAxionEvent->Translate(TVector3(-fCenter.X(), -fCenter.Y(), -fCenter.Z())); + RESTDebug << " ---- " << RESTendl; + RESTDebug << "BoEP: Final Position. X: " << fAxionEvent->GetPosition().X() + << " Y: " << fAxionEvent->GetPosition().Y() << " Z: " << fAxionEvent->GetPosition().Z() + << RESTendl; + RESTDebug << "BoEP: Final Direction. X: " << fAxionEvent->GetDirection().X() + << " Y: " << fAxionEvent->GetDirection().Y() << " Z: " << fAxionEvent->GetDirection().Z() + << RESTendl; + RESTDebug << " ++++ " << RESTendl; +} + +////////////////////////////////////////////////////////////////////////// +/// \brief End of event process. Validate the updated observable number matches total defined observable +/// number +/// +void TRestAxionEventProcess::EndOfEventProcess(TRestEvent* evInput) { + TRestEventProcess::EndOfEventProcess(evInput); + + RESTDebug << "EoEP: Initial Position. X: " << fAxionEvent->GetPosition().X() + << " Y: " << fAxionEvent->GetPosition().Y() << " Z: " << fAxionEvent->GetPosition().Z() + << RESTendl; + RESTDebug << "EoEP: Initial Direction. X: " << fAxionEvent->GetDirection().X() + << " Y: " << fAxionEvent->GetDirection().Y() << " Z: " << fAxionEvent->GetDirection().Z() + << RESTendl; + RESTDebug << " ---- " << RESTendl; + fAxionEvent->Translate(TVector3(fCenter.X(), fCenter.Y(), fCenter.Z())); + fAxionEvent->RotateXY(fCenter, fPitch, fYaw); + RESTDebug << "EoEP: Final Position. X: " << fAxionEvent->GetPosition().X() + << " Y: " << fAxionEvent->GetPosition().Y() << " Z: " << fAxionEvent->GetPosition().Z() + << RESTendl; + RESTDebug << "EoEP: Final Direction. X: " << fAxionEvent->GetDirection().X() + << " Y: " << fAxionEvent->GetDirection().Y() << " Z: " << fAxionEvent->GetDirection().Z() + << RESTendl; + RESTDebug << " ++++ " << RESTendl; +} + +////////////////////////////////////////////////////////////////////////// +/// \brief Pre-defined printer, can be used at the beginning in the +/// implementation of PrintMetadata() +/// +/// Prints process type, name, title, verboselevel, outputlevel, input/output +/// event type, and several separators +void TRestAxionEventProcess::BeginPrintProcess() { + TRestEventProcess::BeginPrintProcess(); + + RESTMetadata << "Center: (" << fCenter.X() << ", " << fCenter.Y() << ", " << fCenter.Z() << ")" + << RESTendl; + RESTMetadata << "Yaw angle (Y-axis): " << fYaw * units("degrees") << " degrees" << RESTendl; + RESTMetadata << "Pitch angle (X-axis): " << fPitch * units("degrees") << " degrees" << RESTendl; + RESTMetadata << " --------------------------- " << RESTendl; + RESTMetadata << " " << RESTendl; +} + +////////////////////////////////////////////////////////////////////////// +/// \brief Adds the footer for PrintMetadata +/// +void TRestAxionEventProcess::EndPrintProcess() { TRestEventProcess::EndPrintProcess(); } diff --git a/src/TRestAxionFieldPropagationProcess.cxx b/src/TRestAxionFieldPropagationProcess.cxx index 45af4012..96b6e4aa 100644 --- a/src/TRestAxionFieldPropagationProcess.cxx +++ b/src/TRestAxionFieldPropagationProcess.cxx @@ -74,11 +74,8 @@ /// #include "TRestAxionFieldPropagationProcess.h" #include -#include "TComplex.h" -#include "TH1F.h" using namespace std; -using namespace REST_Physics; ClassImp(TRestAxionFieldPropagationProcess); @@ -98,56 +95,21 @@ TRestAxionFieldPropagationProcess::TRestAxionFieldPropagationProcess() { Initial /// /// \param cfgFileName A const char* giving the path to an RML file. /// -TRestAxionFieldPropagationProcess::TRestAxionFieldPropagationProcess(char* cfgFileName) { - Initialize(); - - LoadConfig(cfgFileName); -} +TRestAxionFieldPropagationProcess::TRestAxionFieldPropagationProcess(char* cfgFileName) { Initialize(); } /////////////////////////////////////////////// /// \brief Default destructor /// TRestAxionFieldPropagationProcess::~TRestAxionFieldPropagationProcess() { delete fAxionEvent; } -/////////////////////////////////////////////// -/// \brief Function to load the default config in absence of RML input -/// -void TRestAxionFieldPropagationProcess::LoadDefaultConfig() { - SetName(this->ClassName()); - SetTitle("Default config"); -} - -/////////////////////////////////////////////// -/// \brief Function to load the configuration from an external configuration file. -/// -/// If no configuration path is defined in TRestMetadata::SetConfigFilePath -/// the path to the config file must be specified using full path, absolute or relative. -/// -/// \param cfgFileName A const char* giving the path to an RML file. -/// \param name The name of the specific metadata. It will be used to find the -/// corresponding TRestAxionFieldPropagationProcess section inside the RML. -/// -void TRestAxionFieldPropagationProcess::LoadConfig(std::string cfgFilename, std::string name) { - if (LoadConfigFromFile(cfgFilename, name)) LoadDefaultConfig(); -} - /////////////////////////////////////////////// /// \brief Function to initialize input/output event members and define the section name /// -/// 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 TRestAxionFieldPropagationProcess::Initialize() { SetSectionName(this->ClassName()); SetLibraryVersion(LIBRARY_VERSION); - mpfr::mpreal::set_default_prec(mpfr::digits2bits(30)); - fAxionEvent = new TRestAxionEvent(); - - fFinalNormalPlan = TVector3(); - fFinalPositionPlan = TVector3(); - fDistance = 0.0; } /////////////////////////////////////////////// @@ -157,956 +119,36 @@ void TRestAxionFieldPropagationProcess::Initialize() { void TRestAxionFieldPropagationProcess::InitProcess() { RESTDebug << "Entering ... TRestAxionGeneratorProcess::InitProcess" << RESTendl; - fAxionMagneticField = (TRestAxionMagneticField*)this->GetMetadata("TRestAxionMagneticField"); + fField = (TRestAxionMagneticField*)this->GetMetadata("TRestAxionMagneticField"); - if (!fAxionMagneticField) { + if (!fField) { RESTError << "TRestAxionFieldPropagationprocess. Magnetic Field was not defined!" << RESTendl; exit(0); } - - fAxionBufferGas = (TRestAxionBufferGas*)this->GetMetadata("TRestAxionBufferGas"); - - if (!fAxionBufferGas) { - RESTError << "TRestAxionBufferGas. Cannot access the buffer gas" << RESTendl; - exit(0); - } - - fAxionPhotonConversion = new TRestAxionPhotonConversion(); - fAxionPhotonConversion->SetBufferGas(fAxionBufferGas); -} - -/////////////////////////////////////////////// -/// \brief This method will translate the vector with direction `dir` starting at position `pos` to the plane -/// defined by the normal vector plane, `n` that contains the point `a` in the plane. -/// -/// This method has been migrated to TRestPhysics, and it will be accesible through REST_Physics -/* -TVector3 TRestAxionFieldPropagationProcess::MoveToPlane(TVector3 pos, TVector3 dir, TVector3 n, TVector3 a) { - if (n * dir == 0) { - RESTError << "The vector is parallel to the plane!!" << endl; - RESTError << "Position will not be translated" << endl; - } else { - Double_t t = (n * a - n * pos) / (n * dir); - - return pos + t * dir; - } - return pos; -}*/ - -/////////////////////////////////////////////// -/// \brief This method is OBSOLETE. -/// -/// It is more intuitive to define the vector position and direction, and -/// plane vector and point. Then use this information to find the intersection. Instead of defining a -/// component and an impact factor. -/* -TVector3 TRestAxionFieldPropagationProcess::MoveToPlan(TVector3 pos, TVector3 dir, Double_t f, Int_t i) { - if (dir[i] == 0) RESTError << "The component of the direction you chose is equal to 0 " << endl; - - Double_t t = (f - pos[i]) / dir[i]; - pos[i] = f; - pos[(i + 1) % 3] = pos[(i + 1) % 3] + t * dir[(i + 1) % 3]; - pos[(i + 2) % 3] = pos[(i + 2) % 3] + t * dir[(i + 2) % 3]; - - return pos; -} -*/ - -/////////////////////////////////////////////// -/// \brief This method is OBSOLETE. -/// -/// Re-implementation of MoveToPlane -/* -TVector3 TRestAxionFieldPropagationProcess::FinalPositionInPlan(TVector3 pos, TVector3 dir, - TVector3 normalPlan, TVector3 pointPlan) { - if (normalPlan.Dot(dir) == 0) return pos; - - Double_t t, d; - d = -normalPlan.Dot(pointPlan); - t = -(normalPlan.Dot(pos) + d) / (normalPlan.Dot(dir)); - - pos = pos + t * dir; - - return pos; -} -*/ - -/////////////////////////////////////////////// -/// \brief This method is OBSOLETE -/// -/// Re-implementation in MoveByDistance. Just because the method name is more intuitive using BYDISTANCE -/* -TVector3 TRestAxionFieldPropagationProcess::MoveToFinalDistance(TVector3 pos, TVector3 dir, - Double_t distance) { - Double_t t = distance / dir.Mag(); - pos = pos + t * dir; - - return pos; -} -*/ - -/////////////////////////////////////////////// -/// \brief This method transports a position `pos` by a distance `d` in the direction defined by `dir`. -/// -/// This method has been migrated to TRestPhysics, and it will be accesible through REST_Physics -/* -TVector3 TRestAxionFieldPropagationProcess::MoveByDistance(TVector3 pos, TVector3 dir, Double_t d) { - return pos + d * dir.Unit(); -} */ - -bool TRestAxionFieldPropagationProcess::IsInBoundedPlan(TVector3 pos, Int_t i, Int_t p) { - Double_t minCond1, maxCond1, minCond2, maxCond2; - Int_t j, k; - - /* -if (i == 0) { - minCond1 = (fAxionMagneticField->GetYmin())[p]; - maxCond1 = (fAxionMagneticField->GetYmax())[p]; - minCond2 = (fAxionMagneticField->GetZmin())[p]; - maxCond2 = (fAxionMagneticField->GetZmax())[p]; - j = 1; - k = 2; -} -if (i == 1) { - minCond1 = (fAxionMagneticField->GetXmin())[p]; - maxCond1 = (fAxionMagneticField->GetXmax())[p]; - minCond2 = (fAxionMagneticField->GetZmin())[p]; - maxCond2 = (fAxionMagneticField->GetZmax())[p]; - j = 0; - k = 2; -} -if (i == 2) { - minCond1 = (fAxionMagneticField->GetXmin())[p]; - maxCond1 = (fAxionMagneticField->GetXmax())[p]; - minCond2 = (fAxionMagneticField->GetYmin())[p]; - maxCond2 = (fAxionMagneticField->GetYmax())[p]; - j = 0; - k = 1; -} - - bool cond1 = (pos[j] <= maxCond1 && pos[j] >= minCond1); - bool cond2 = (pos[k] <= maxCond2 && pos[k] >= minCond2); - - return cond1 && cond2; - */ - return false; -} - -std::vector TRestAxionFieldPropagationProcess::InOut(std::vector bounds, TVector3 dir) { - Int_t i = 0; - - TVector3 in = bounds[0]; - TVector3 out = bounds[1]; - - while (i < 3 && dir[i] * (out[i] - in[i]) >= 0) i = i + 1; - - if (i < 3) { - bounds[0] = out; - bounds[1] = in; - } - - return bounds; -} - -/////////////////////////////////////////////// -/// \brief Finds the in/out particle trajectory boundaries for a particular magnetic volume. -/// -/// This method checks if the particle (with the initial position `pos` and direction `dir`) passes through -/// the magnetic field region specified by the input parameter p. It is done by searching for the points where -/// the particle trajectory intersects the boundary planes of that region. If two such points (entry point and -/// exit point) are found, their coordinates are stored in the vector boundaries. In the example shown in Fig. -/// 1 these points are: IN 1 and OUT 1 for the region #1 and IN2 and OUT 2 for the region #2. -/* This method has been moved to TRestAxionMagneticField::GetVolumeBoundaries -std::vector TRestAxionFieldPropagationProcess::FindBoundariesOneVolume(TVector3 pos, TVector3 dir, - Int_t p) { - std::vector boundaries; -TVector3 in; -TVector3 out; - -Int_t i = 0; -Int_t j = 0; - -TVectorD f(6); -f[0] = (fAxionMagneticField->GetXmin())[p]; -f[1] = (fAxionMagneticField->GetXmax())[p]; -f[2] = (fAxionMagneticField->GetYmin())[p]; -f[3] = (fAxionMagneticField->GetYmax())[p]; -f[4] = (fAxionMagneticField->GetZmin())[p]; -f[5] = (fAxionMagneticField->GetZmax())[p]; - -Int_t nFace = 0; - -while (nFace < 6 && i <= 0) { - if (dir[Int_t(nFace / 2)] != 0) { - if ((pos[Int_t(nFace / 2)] - f[nFace]) * dir[Int_t(nFace / 2)] <= 0) { - in = MoveToPlan(pos, dir, f[nFace], Int_t(nFace / 2)); - if (IsInBoundedPlan(in, Int_t(nFace / 2), p)) i = i + 1; - } - } - nFace = nFace + 1; -} - -if (i == 1) { - while (nFace < 6 && j <= 0) { - if (dir[Int_t(nFace / 2)] != 0) { - if ((pos[Int_t(nFace / 2)] - f[nFace]) * dir[Int_t(nFace / 2)] <= 0) { - out = MoveToPlan(pos, dir, f[nFace], Int_t(nFace / 2)); - if (IsInBoundedPlan(out, Int_t(nFace / 2), p)) j = j + 1; - } - } - nFace = nFace + 1; - } -} - -if (i + j == 2 && in != out) { - boundaries.push_back(in); - boundaries.push_back(out); - return boundaries; -} - -else - return boundaries; - return boundaries; -} -*/ - -/* This method has been moved to TRestAxionMagneticField::GetFieldBoundaries -std::vector TRestAxionFieldPropagationProcess::FieldBoundary(std::vector boundaries, - Double_t minStep) { - std::vector boundariesField; - - TVector3 in = boundaries[0]; - TVector3 out = boundaries[1]; - TVector3 diff = out - in; - - Int_t N = 10; - - Int_t i = 0; - - while (1.0 / Double_t(N) > minStep) { - while ((fAxionMagneticField->GetMagneticField(in[0], in[1], in[2])) == TVector3(0, 0, 0) && i < N) { - in = in + 1.0 / Double_t(N) * diff; - i = i + 1; - } - - if (i == N) return boundariesField; - - in = in - 1.0 / Double_t(N) * diff; - N = 10 * N; - i = 0; - } - - in = in + 10.0 / Double_t(N) * diff; - boundariesField.push_back(in); - - N = 10; - i = 0; - - while (1.0 / Double_t(N) > minStep) { - while ((fAxionMagneticField->GetMagneticField(out[0], out[1], out[2]) == TVector3(0, 0, 0)) && - i < N) { - out = out - 1.0 / Double_t(N) * diff; - i = i + 1; - } - - out = out + 1.0 / Double_t(N) * diff; - N = 10 * N; - i = 0; - } - - out = out - 10.0 / Double_t(N) * diff; - boundariesField.push_back(out); - - return boundariesField; -} -*/ - -// This method is obsolete because it uses methods FindBoundariesOneVolume and FieldBoundary that are obsolete -// and replaced by the methods TRestAxionMagneticField::GetVolumeBoundaries and -// TRestAxionMagneticField::GetFieldBoundaries. -// The NEW version of this method is given below and uses TRestAxionMagneticField::GetFieldBoundaries -/* -std::vector> TRestAxionFieldPropagationProcess::FindFieldBoundaries(Double_t minStep) { - std::vector> boundaryCollection; - std::vector> boundaryFinalCollection; - -if (minStep == -1) minStep = 0.01; -std::vector buffVect; -TVector3 boundaryIn; -TVector3 boundaryOut; - -TVector3 posInitial = *(fInputAxionEvent->GetPosition()); -TVector3 direction = *(fInputAxionEvent->GetDirection()); -direction = direction.Unit(); - -if (direction == TVector3(0, 0, 0)) // No moves - return boundaryCollection; - -if ((fAxionMagneticField->GetXmin()).size() == 0) - RESTError << " The magnetic field has not been loaded " << endl; - -// Find global boundaries volume - -Int_t N = (fAxionMagneticField->GetXmin()).size(); - -for (Int_t p = 0; p < N; p++) { - buffVect = FindBoundariesOneVolume(posInitial, direction, p); - if (buffVect.size() == 2) { - buffVect = InOut(buffVect, direction); - boundaryCollection.push_back(buffVect); - buffVect.clear(); - } -} - -debug << "+------------------------+" << endl; -debug << " Number of volume boundaries : " << 2 * boundaryCollection.size() << endl; -debug << "+------------------------+" << endl; - -debug << "+------------------------+" << endl; -for (Int_t p = 0; p < boundaryCollection.size(); p++) { - debug << "for" << p << " in : (" << boundaryCollection[p][0].X() << "," - << boundaryCollection[p][0].Y() << "," << boundaryCollection[p][0].Z() << ")" << endl; - debug << "for" << p << " out : (" << boundaryCollection[p][1].X() << "," - << boundaryCollection[p][1].Y() << "," << boundaryCollection[p][1].Z() << ")" << endl; -} -debug << "+------------------------+" << endl; - -// Find precise boundaries field - -for (Int_t i = 0; i < boundaryCollection.size(); i++) { - buffVect = FieldBoundary(boundaryCollection[i], minStep); - if (buffVect.size() == 2) { - boundaryFinalCollection.push_back(buffVect); - buffVect.clear(); - } -} - -boundaryCollection.clear(); - -debug << "+------------------------+" << endl; -debug << " Number of field boundaries : " << 2 * boundaryFinalCollection.size() << endl; -debug << "+------------------------+" << endl; - -debug << "+------------------------+" << endl; -for (Int_t p = 0; p < boundaryFinalCollection.size(); p++) { - debug << "for" << p << " in : (" << boundaryFinalCollection[p][0].X() << "," - << boundaryFinalCollection[p][0].Y() << "," << boundaryFinalCollection[p][0].Z() << ")" << endl; - debug << "for" << p << " out : (" << boundaryFinalCollection[p][1].X() << "," - << boundaryFinalCollection[p][1].Y() << "," << boundaryFinalCollection[p][1].Z() << ")" << endl; -} -debug << "+------------------------+" << endl; - - return boundaryFinalCollection; -} -*/ - -// This is a NEW version of FindFieldBoundaries method and replaces the obsolete version of this method given -// above -/////////////////////////////////////////////// -/// \brief Finds the boundary points for each segment along the particle trajectory where the **transversal** -/// component -/// of the magnetic field in not zero. -/// -/// This method goes over all magnetic field regions and checks for each region if the particle (with the -/// initial -/// position `posInitial` and direction `direction`) passes through this region. If yes, it determines the -/// boundary -/// points of the trajectory segment through this region between which the **transversal** component of the -/// magnetic -/// field is not zero. The method returns the vector where each element contains a pair of boundary points, -/// i.e., the starting and ending points of one segment of the particle trajectory along which the -/// **transversal** -/// component of the magnetic field is not zero. -/// -std::vector> TRestAxionFieldPropagationProcess::FindFieldBoundaries(Double_t minStep) { - std::vector> boundaryFinalCollection; - - if (minStep == -1) minStep = 0.01; - std::vector buffVect; - - TVector3 posInitial = *(fAxionEvent->GetPosition()); - TVector3 direction = *(fAxionEvent->GetDirection()); - direction = direction.Unit(); - - if (direction == TVector3(0, 0, 0)) // No moves - return boundaryFinalCollection; - - if (fAxionMagneticField->GetNumberOfVolumes() == 0) - RESTError << " The magnetic field has not been loaded " << RESTendl; - - // Find field boundaries volume - - Int_t N = fAxionMagneticField->GetNumberOfVolumes(); - - for (Int_t p = 0; p < N; p++) { - buffVect.clear(); - buffVect = fAxionMagneticField->GetFieldBoundaries(posInitial, direction, p); - if (buffVect.size() == 2) { - boundaryFinalCollection.push_back(buffVect); - } - } - - RESTDebug << "+------------------------+" << RESTendl; - RESTDebug << " Number of field boundaries : " << 2 * boundaryFinalCollection.size() << RESTendl; - RESTDebug << "+------------------------+" << RESTendl; - - RESTDebug << "+------------------------+" << RESTendl; - for (Int_t p = 0; p < boundaryFinalCollection.size(); p++) { - RESTDebug << "for volume " << p << " in : (" << boundaryFinalCollection[p][0].X() << "," - << boundaryFinalCollection[p][0].Y() << "," << boundaryFinalCollection[p][0].Z() << ")" - << RESTendl; - RESTDebug << "for volume " << p << " out : (" << boundaryFinalCollection[p][1].X() << "," - << boundaryFinalCollection[p][1].Y() << "," << boundaryFinalCollection[p][1].Z() << ")" - << RESTendl; - } - RESTDebug << "+------------------------+" << RESTendl; - - return boundaryFinalCollection; -} - -/* This method is now OBSOLETE. It seems to me there is a problem here, there are two directional vectors -defined. -The one defined by in/out coordinates, and the one defined by axion direction. - -This method will be substituted by -std::vector TRestAxionMagneticField::GetTransversalComponentAlongPath(TVector3 from, TVector3 to, -Double_t dl, Int_t Nmax ); - -TVectorD TRestAxionFieldPropagationProcess::GetFieldVector(TVector3 in, TVector3 out, Int_t N) { - if (N == 0) N = TMath::Power(10, 4); - - TVectorD Bt(N); - TVector3 B; - TVector3 direction = *(fInputAxionEvent->GetDirection()); - - TVector3 differential = out - in; - - B = fAxionMagneticField->GetMagneticField(in[0], in[1], in[2]); - Bt[0] = abs(B.Perp(direction)); - - for (Int_t i = 1; i < N; i++) { - in = in + differential * (1.0 / Double_t(N - 1)); - B = fAxionMagneticField->GetMagneticField(in[0], in[1], in[2]); - Bt[i] = abs(B.Perp(direction)); - } - - return Bt; -} */ - -/// \brief Prints variables of the ComplexReal type, i.e, complex numbers -/// -void TRestAxionFieldPropagationProcess::PrintComplex(ComplexReal p) { - RESTDebug << p.real << " + " << p.img << "i" << RESTendl; -} - -/// \brief Calculates amplitudes of the axion field, parallel component of the photon field and orthogonal -/// component of the photon field (parallel and orthogonal are with respect to the transverse component of the -/// magnetic field) -/// after passing one segment of the particle trajectory along which it is assumed -/// that the transverse component of the magnetic field with respect to the particle propagation direction -/// is not equal zero. The segment is divided into a series of subsegments. It is assumed that the magnitude -/// and direction of the -/// transverse component of the magnetic field are constant along each subsegment, but they can change from -/// one subsegment -/// to another. The calculations are based on the procedure described in the Section 4 of the internal report -/// "Axion-photon conversion in the external magnetic fields" written by B. Lakic and K. Jakovcic. -/// For each subsegment, the values of certain parameters (e.g. theta, lambda) are calculated first and then -/// the method -/// `CalculateAmplitudesInSubsegment` is called to calculate the amplitudes after passing that subsegment, -/// i.e., to calculate the amplitudes -/// at the end of that subsegment. This procedure is repeated for each subsegment until -/// the end of segment is reached. The values of the amplitudes at the end of one subsegment are used to -/// calculate the initial -/// values of these amplitudes for the next subsegment by using equations (4.4)-(4.6). -/// -/// NOTE: The amplitudes are calculated for the axion-photon coupling constant g_agg = 10^-10 GeV-1 -/// The length of the subsegment is defined by variable `step`. It is currently set to 200 mm. - -void TRestAxionFieldPropagationProcess::CalculateAmplitudesInSegment( - ComplexReal& faxionAmplitude, ComplexReal& fparallelPhotonAmplitude, - ComplexReal& forthogonalPhotonAmplitude, TVector3& averageBT, mpfr::mpreal axionMass, - mpfr::mpreal photonMass, mpfr::mpreal Ea, TVector3 from, TVector3 to, mpfr::mpreal CommonPhase, - mpfr::mpreal OrthogonalPhase) { - mpfr::mpreal g_agg = 1.0e-10; // axion-photon coupling constant in GeV-1 - mpfr::mpreal TeslaineV = 195.3; // conversion factor from Tesla to eV^2 - Double_t segment_length = (to - from).Mag(); // the length of the entire segment - TVector3 subsegment_start, - subsegment_end; // coordinates of the starting and ending points of one subsegment - Double_t subsegment_length; // the length of one subsegment - ComplexReal subsegment_A0_par, subsegment_A0_ort; // values of the parallel and orthogonal photon - // amplitudes at the end of one subsegment (that is - // also the beginning of the next subsegment - TVector3 - averageBT_0; // transverse component of the average magnetic field vector in the previous subsegment - - Double_t step = 200.0; // the length of the subsegment - Double_t BTmag; // average magnitude of the transverse component of the magnetic field in one subsegment - // (given in Tesla) - Double_t BTangle; // angle between the transverse component of the average magnetic field in one - // subsegment and the one in the previous subsegment - - subsegment_start = from; - subsegment_A0_par = fparallelPhotonAmplitude; // parallel photon amplitude at the beginning of the - // segment, i.e., at the point `from` - subsegment_A0_ort = forthogonalPhotonAmplitude; // orthogonal photon amplitude at the beginning of the - // segment, i.e., at the point `from` - averageBT_0 = averageBT; - TVector3 dir = (to - from).Unit(); // direction of the particle propagation - - while ((subsegment_start - from).Mag() < segment_length) { // loop over subsegments in one segment - subsegment_end = subsegment_start + step * dir; - if ((subsegment_end - from).Mag() >= segment_length) subsegment_end = to; - subsegment_length = (subsegment_end - subsegment_start).Mag(); - subsegment_length = subsegment_length / 1000.0; // default REST units are mm - - // calculation of the average magnitude of the transverse magnetic field along the subsegment, i.e., - // between coordinates `subsegment_start` and `subsegment_end` - BTmag = - fAxionMagneticField->GetTransversalFieldAverage(subsegment_start, subsegment_end); // in Tesla - - // calculation of the transverse component of the average magnetic field vector along the subsegment, - // i.e., between coordinates `subsegment_start` and `subsegment_end` - averageBT = fAxionMagneticField->GetFieldAverageTransverseVector(subsegment_start, subsegment_end); - - // calculation of the angle between the transverse component of the average magnetic field in one - // subsegment and the one in the previous subsegment - if ((averageBT_0.Mag() == 0.0) || (averageBT.Mag() == 0.0)) - BTangle = 0.0; - else - BTangle = averageBT.Angle(averageBT_0); - if (averageBT.Mag() != 0.0) averageBT_0 = averageBT; - - // calculation of the initial values of the parallel and orthogonal photon amplitudes at the beginning - // of the subsegment (at the point `subsegment_start`) - // from the values at the end of the previous subsegment (also at the point `subsegment_start`) - fparallelPhotonAmplitude = ComplexAddition(ComplexProduct(cos(BTangle), subsegment_A0_par), - ComplexProduct(sin(BTangle), subsegment_A0_ort)); - forthogonalPhotonAmplitude = ComplexAddition(ComplexProduct(-sin(BTangle), subsegment_A0_par), - ComplexProduct(cos(BTangle), subsegment_A0_ort)); - - // calculation of the parameters theta and lambda for the subsegment - mpfr::mpreal term_1 = 2 * (Ea * 1000.0) * (g_agg * 1.0e-9) * (BTmag * TeslaineV); // in eV^2 - mpfr::mpreal term_2 = axionMass * axionMass - photonMass * photonMass; // in ev^2 - mpfr::mpreal theta = 0.5 * atan(term_1 / term_2); - mpfr::mpreal lambda = sqrt(term_1 * term_1 + term_2 * term_2) / (4. * Ea * 1000.0); // in eV - - RESTDebug << "+--------------------------------------------------------------------------+" - << RESTendl; - RESTDebug << " CalculateAmplitudesInSegment method: Parameter summary" << RESTendl; - RESTDebug << RESTendl << "segment length = " << segment_length << " mm" << RESTendl; - RESTDebug << RESTendl << "subsegment_start: (" << subsegment_start.x() << ", " << subsegment_start.y() - << ", " << subsegment_start.z() << ") " - << " mm" << RESTendl; - RESTDebug << "subsegment_end: ( " << subsegment_end.x() << ", " << subsegment_end.y() << ", " - << subsegment_end.z() << ") " - << " mm" << RESTendl; - RESTDebug << "subsegment length = " << subsegment_length << " m" << RESTendl << RESTendl; - - RESTDebug - << " average magnitude of the transverse component of the magnetic field in the subsegment : " - << BTmag << " T" << RESTendl; - RESTDebug - << " angle of the transverse component of the average magnetic field in the subsegment with " - "respect to the previous subsegment : " - << BTangle << " rad" << RESTendl; - RESTDebug << " g_agg : " << g_agg << " GeV-1" << RESTendl; - RESTDebug << " Theta : " << theta << RESTendl; - RESTDebug << " lambda : " << lambda << " eV" << RESTendl; - RESTDebug << " subsegment_A0_par : "; - PrintComplex(subsegment_A0_par); - RESTDebug << " subsegment_A0_ort : "; - PrintComplex(subsegment_A0_ort); - RESTDebug << " BEFORE calculating in subsegment: paralell photon component amplitude : "; - PrintComplex(fparallelPhotonAmplitude); - RESTDebug << " BEFORE calculating in subsegment: orthogonal photon component amplitude : "; - PrintComplex(forthogonalPhotonAmplitude); - RESTDebug << "+--------------------------------------------------------------------------+" - << RESTendl; - - CalculateAmplitudesInSubsegment(faxionAmplitude, fparallelPhotonAmplitude, forthogonalPhotonAmplitude, - theta, lambda, subsegment_length, CommonPhase, OrthogonalPhase); - subsegment_A0_par = fparallelPhotonAmplitude; // parallel photon amplitude at the end of the - // subsegment, i.e. at the point `subsegment_end` - subsegment_A0_ort = forthogonalPhotonAmplitude; // orthogonal photon amplitude at the end of the - // subsegment, i.e. at the point `subsegment_end` - RESTDebug << RESTendl << " AFTER calculating in subsegment: subsegment_A0_par : "; - PrintComplex(subsegment_A0_par); - RESTDebug << " AFTER calculating in subsegment: subsegment_A0_ort : "; - PrintComplex(subsegment_A0_ort); - RESTDebug << " AFTER calculating in subsegment: paralell photon component amplitude : "; - PrintComplex(fparallelPhotonAmplitude); - RESTDebug << " AFTER calculating in subsegment: orthogonal photon component amplitude : "; - PrintComplex(forthogonalPhotonAmplitude); - subsegment_start = subsegment_end; - } -} - -/// \brief Calculates amplitudes of the axion field, parallel component of the photon field and orthogonal -/// component -/// of the photon field after passing one subsegment of the particle trajectory along which it is assumed -/// that the transverse component of the magnetic field is constant. It uses equations (3.15) and (3.38) given -/// in the internal report "Axion-photon conversion in the external magnetic fields" written by B. Lakic and -/// K. Jakovcic. -/// Also, amplitudes are of ComplexReal type, which stores complex numbers based on mpreal wrapper to allow -/// precise -/// calculation of small values. At present, the precision is set to 30 digits, so that we can still -/// calculate -/// numbers such as : 1.0 - 1.e-30 -/// NOTE: The amplitudes are calculated for the axion-photon coupling constant g_agg = 10^-10 GeV-1 -void TRestAxionFieldPropagationProcess::CalculateAmplitudesInSubsegment( - ComplexReal& faxionAmplitude, ComplexReal& fparallelPhotonAmplitude, - ComplexReal& forthogonalPhotonAmplitude, mpfr::mpreal theta, mpfr::mpreal lambda, Double_t length, - mpfr::mpreal CommonPhase, mpfr::mpreal OrthogonalPhase) { - // lambda is given in eV, theta has no dimension, length is in m, Common phase and Orthogonal phase are in - // eV - - mpfr::mpreal::set_default_prec(mpfr::digits2bits(30)); - cout.precision(30); - // setting initial parameters - ComplexReal a0 = - faxionAmplitude; // axion field amplitude initial value at the beginning of the subsegment - ComplexReal A0_par = fparallelPhotonAmplitude; // photon field parallel component amplitude initial value - // at the beginning of the subsegment - ComplexReal A0_ort = forthogonalPhotonAmplitude; // photon field orthogonal component amplitude initial - // value at the beginning of the subsegment - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " CalculateAmplitudesInSubsegment method: Parameter summary" << RESTendl; - RESTDebug << " Theta : " << theta << RESTendl; - RESTDebug << " lambda : " << lambda << " eV" << RESTendl; - RESTDebug << " subsegment length : " << length << " m" << RESTendl; - RESTDebug << " Common phase : " << CommonPhase << " eV" << RESTendl; - RESTDebug << " orthogonal photon component phase : " << OrthogonalPhase << " eV" << RESTendl; - RESTDebug << " a0 : "; - PrintComplex(a0); - RESTDebug << " A0_par : "; - PrintComplex(A0_par); - RESTDebug << " A0_ort : "; - PrintComplex(A0_ort); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - - // setting auxillary parameters used in calculations - mpfr::mpreal cos2_theta = cos(theta) * cos(theta); - mpfr::mpreal sin2_theta = sin(theta) * sin(theta); - mpfr::mpreal sin_2theta = sin(2.0 * theta); - - mpfr::mpreal l = length * PhMeterIneV; // length in eV-1 - - mpfr::mpreal phi = lambda * l; - ComplexReal exp_lambdaPlusZ = SetComplexReal(cos(-phi), sin(-phi)); - ComplexReal exp_lambdaMinusZ = SetComplexReal(cos(phi), sin(phi)); - - mpfr::mpreal phi1 = CommonPhase * l; - ComplexReal exp_CommonPhase = SetComplexReal(cos(phi1), sin(phi1)); - - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " Intermediate calculations" << RESTendl; - RESTDebug << " cos^(2)_theta : " << cos2_theta << RESTendl; - RESTDebug << " sin^(2)_theta : " << sin2_theta << RESTendl; - RESTDebug << " sin(2*theta) : " << sin_2theta << RESTendl; - RESTDebug << " subsegment length : " << length << " m" << RESTendl; - RESTDebug << " l : " << l << " eV-1" << RESTendl; - RESTDebug << " phi : " << phi << RESTendl; - RESTDebug << " exp_lambdaPlusZ : "; - PrintComplex(exp_lambdaPlusZ); - RESTDebug << " exp_lambdaMinusZ : "; - PrintComplex(exp_lambdaMinusZ); - RESTDebug << " phi1 : " << phi1 << RESTendl; - RESTDebug << " exp_CommonPhase : "; - PrintComplex(exp_CommonPhase); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - - // calculation of the photon parallel component amplitude - ComplexReal A_term_1_1 = ComplexProduct(cos2_theta, exp_lambdaPlusZ); - ComplexReal A_term_1_2 = ComplexProduct(sin2_theta, exp_lambdaMinusZ); - ComplexReal A_sum_1 = ComplexAddition(A_term_1_1, A_term_1_2); - ComplexReal A_term_1 = ComplexProduct(A0_par, A_sum_1); - ComplexReal A_sum_2 = ComplexSubstraction(exp_lambdaPlusZ, exp_lambdaMinusZ); - ComplexReal temp = ComplexProduct(sin_2theta * 0.5, a0); - ComplexReal A_term_2 = ComplexProduct(temp, A_sum_2); - ComplexReal A_sum = ComplexAddition(A_term_1, A_term_2); - fparallelPhotonAmplitude = ComplexProduct(exp_CommonPhase, A_sum); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " Intermediate calculations for the photon parallel component amplitude" << RESTendl; - RESTDebug << " A_term_1_1 : "; - PrintComplex(A_term_1_1); - RESTDebug << " A_term_1_2 : "; - PrintComplex(A_term_1_2); - RESTDebug << " A_sum_1 : "; - PrintComplex(A_sum_1); - RESTDebug << " A_term_1 : "; - PrintComplex(A_term_1); - RESTDebug << " A_sum_2 : "; - PrintComplex(A_sum_2); - RESTDebug << " A_term_2 : "; - PrintComplex(A_term_2); - RESTDebug << " A_sum : "; - PrintComplex(A_sum); - RESTDebug << " parallelPhotonAmplitude : "; - PrintComplex(fparallelPhotonAmplitude); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - - // calculation of the axion amplitude - ComplexReal a_sum_1 = A_sum_2; - temp = ComplexProduct(sin_2theta * 0.5, A0_par); - ComplexReal a_term_1 = ComplexProduct(temp, a_sum_1); - ComplexReal a_term_2_1 = ComplexProduct(sin2_theta, exp_lambdaPlusZ); - ComplexReal a_term_2_2 = ComplexProduct(cos2_theta, exp_lambdaMinusZ); - ComplexReal a_sum_2 = ComplexAddition(a_term_2_1, a_term_2_2); - ComplexReal a_term_2 = ComplexProduct(a0, a_sum_2); - ComplexReal a_sum = ComplexAddition(a_term_1, a_term_2); - faxionAmplitude = ComplexProduct(exp_CommonPhase, a_sum); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " Intermediate calculations for the axion amplitude" << RESTendl; - RESTDebug << " a_sum_1 : "; - PrintComplex(a_sum_1); - RESTDebug << " a_term_1 : "; - PrintComplex(a_term_1); - RESTDebug << " a_term_2_1 : "; - PrintComplex(a_term_2_1); - RESTDebug << " a_term_2_2 : "; - PrintComplex(a_term_2_2); - RESTDebug << " a_sum_2 : "; - PrintComplex(a_sum_2); - RESTDebug << " a_term_2 : "; - PrintComplex(a_term_2); - RESTDebug << " a_sum : "; - PrintComplex(a_sum); - RESTDebug << " axionAmplitude : "; - PrintComplex(faxionAmplitude); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - - // calculation of the photon orthogonal component amplitude - mpfr::mpreal phi2 = OrthogonalPhase * l; - ComplexReal exp_OrthogonalPhase = SetComplexReal(cos(phi2), sin(phi2)); - forthogonalPhotonAmplitude = ComplexProduct(A0_ort, exp_OrthogonalPhase); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " Intermediate calculations for the photon orthogonal component amplitude" << RESTendl; - RESTDebug << " phi2 : " << phi2 << RESTendl; - RESTDebug << " exp_OrthogonalPhase : "; - PrintComplex(exp_OrthogonalPhase); - RESTDebug << " orthogonalPhotonAmplitude : "; - PrintComplex(forthogonalPhotonAmplitude); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; -} - -/// \brief Calculates amplitudes of the axion field, parallel component of the photon field and orthogonal -/// component -/// of the photon field after passing one segment of the particle trajectory along which the transverse -/// component -/// of the magnetic field is ZERO. It uses equation (3.15) given in the internal report "Axion-photon -/// conversion -/// in the external magnetic fields" written by B. Lakic and K. Jakovcic. -/// Also, amplitudes are of ComplexReal type, which stores complex numbers based on mpreal wrapper to allow -/// precise -/// calculation of small values. At present, the precision is set to 30 digits, so that we can still -/// calculate -/// numbers such as : 1.0 - 1.e-30 -void TRestAxionFieldPropagationProcess::PropagateWithoutBField(ComplexReal& faxionAmplitude, - ComplexReal& fparallelPhotonAmplitude, - ComplexReal& forthogonalPhotonAmplitude, - mpfr::mpreal axionMass, - mpfr::mpreal photonMass, mpfr::mpreal Ea, - TVector3 from, TVector3 to) { - mpfr::mpreal::set_default_prec(mpfr::digits2bits(30)); - cout.precision(30); - mpfr::mpreal axionPhase = Ea * 1000.0 - (axionMass * axionMass) / (2. * Ea * 1000.0); // in eV - mpfr::mpreal photonPhase = Ea * 1000.0 - (photonMass * photonMass) / (2. * Ea * 1000.0); // in eV - Double_t length = (to - from).Mag(); - length = length / 1000.0; // default REST units are mm - mpfr::mpreal l = length * PhMeterIneV; // length in eV-1 - - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " Propagation without B field: " << RESTendl; - RESTDebug << " INITIAL VALUES " << RESTendl; - RESTDebug << " axionAmplitude : "; - PrintComplex(faxionAmplitude); - RESTDebug << " parallelPhotonAmplitude : "; - PrintComplex(fparallelPhotonAmplitude); - RESTDebug << " orthogonalPhotonAmplitude : "; - PrintComplex(forthogonalPhotonAmplitude); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - - mpfr::mpreal axionphi = axionPhase * l; - ComplexReal exp_axionPhase = SetComplexReal(cos(axionphi), sin(axionphi)); - faxionAmplitude = ComplexProduct(faxionAmplitude, exp_axionPhase); - - mpfr::mpreal photonphi = photonPhase * l; - ComplexReal exp_photonPhase = SetComplexReal(cos(photonphi), sin(photonphi)); - fparallelPhotonAmplitude = ComplexProduct(fparallelPhotonAmplitude, exp_photonPhase); - forthogonalPhotonAmplitude = ComplexProduct(forthogonalPhotonAmplitude, exp_photonPhase); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " Intermediate calculations" << RESTendl; - RESTDebug << " axionPhase : " << axionPhase << RESTendl; - RESTDebug << " axionphi : " << axionphi << RESTendl; - RESTDebug << " exp_axionPhase : "; - PrintComplex(exp_axionPhase); - RESTDebug << "Norm2(exp_axionPhase) = " << Norm2(exp_axionPhase) << RESTendl; - RESTDebug << " photonPhase : " << photonPhase << RESTendl; - RESTDebug << " photonphi : " << photonphi << RESTendl; - RESTDebug << " exp_photonPhase : "; - PrintComplex(exp_photonPhase); - RESTDebug << "Norm2(exp_photonPhase) = " << Norm2(exp_photonPhase) << RESTendl; - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - RESTDebug << " Propagation without B field: " << RESTendl; - RESTDebug << " FINAL VALUES " << RESTendl; - RESTDebug << " axionAmplitude : "; - PrintComplex(faxionAmplitude); - RESTDebug << " parallelPhotonAmplitude : "; - PrintComplex(fparallelPhotonAmplitude); - RESTDebug << " orthogonalPhotonAmplitude : "; - PrintComplex(forthogonalPhotonAmplitude); - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; } TRestEvent* TRestAxionFieldPropagationProcess::ProcessEvent(TRestEvent* evInput) { - mpfr::mpreal::set_default_prec(mpfr::digits2bits(30)); - cout.precision(30); + // Already done by TRestAxionEventProcess fAxionEvent = (TRestAxionEvent*)evInput; + RESTDebug << "TRestAxionFieldPropagationProcess::ProcessEvent : " << fAxionEvent->GetID() << RESTendl; - TVector3 position = *(fAxionEvent->GetPosition()); - TVector3 direction = *(fAxionEvent->GetDirection()); - Double_t Ea = fAxionEvent->GetEnergy(); - Double_t ma = fAxionEvent->GetMass(); - - faxionAmplitude = SetComplexReal(1.0, 0.0); - fparallelPhotonAmplitude = SetComplexReal(0.0, 0.0); - forthogonalPhotonAmplitude = SetComplexReal(0.0, 0.0); - TVector3 averageBT = TVector3(0.0, 0.0, 0.0); // initial value of the transverse component of the average - // magnetic field before entering the magnetic field region - - RESTDebug << "+------------------------+" << RESTendl; - RESTDebug << "Initial position of the axion input : " << RESTendl; - RESTDebug << "(" << position.X() << "," << position.Y() << "," << position.Z() << ")" << RESTendl; - RESTDebug << "Direction of the axion input : " << RESTendl; - RESTDebug << "(" << direction.X() << "," << direction.Y() << "," << direction.Z() << ")" << RESTendl; - RESTDebug << "Axion energy : " << Ea << RESTendl; - RESTDebug << "Axion mass : " << ma << RESTendl; - RESTDebug << "axion amplitude = " << faxionAmplitude.real << " + " << faxionAmplitude.img << "i" - << RESTendl; - RESTDebug << "parallel photon amplitude = " << fparallelPhotonAmplitude.real << " + " - << fparallelPhotonAmplitude.img << "i" << RESTendl; - RESTDebug << "orthogonal photon amplitude = " << forthogonalPhotonAmplitude.real << " + " - << forthogonalPhotonAmplitude.img << "i" << RESTendl; - RESTDebug << "+------------------------+" << RESTendl; - - std::vector> boundaries; - boundaries = FindFieldBoundaries(); - Int_t NofVolumes = boundaries.size(); - - RESTDebug << "+------------------------+" << RESTendl; - RESTDebug << "Number of magnetic field regions through which the axion passes : " << NofVolumes - << RESTendl; - RESTDebug << "+------------------------+" << RESTendl; - - Double_t probability = 0.; // initial value of the axion to photon conversion probability - mpfr::mpreal axionMass = ma; // in eV - mpfr::mpreal photonMass = 0.0; + std::vector trackBounds = + fField->GetFieldBoundaries(fAxionEvent->GetPosition(), fAxionEvent->GetDirection()); - for (Int_t i = 0; i < NofVolumes; i++) { // loop over segments of trajectory where B is not equal to zero - Int_t id = fAxionMagneticField->GetVolumeIndex(boundaries[i][0]); - if (id < 0) { - RESTWarning << "TRestAxionFieldPropagationProcess::ProcessEvent position is outside any volume. " - "Setting photon mass to 0." - << RESTendl; - photonMass = 0.0; // in eV - } else { - Double_t mphoton = fAxionMagneticField->GetPhotonMass( - id, - Ea); // It returns the effective photon mass in eV at the corresponding magnetic volume id. - photonMass = mphoton; - } - RESTDebug << "Volume ID = " << id << " photon mass = " << photonMass << RESTendl; - RESTDebug << "Volume ID = " << id << " axion mass = " << axionMass << RESTendl; - RESTDebug << "Volume ID = " << id << " axion energy = " << Ea << RESTendl; - - // calculating common phase for the axion field and parallel component of the photon field (eqs. (4.1) - // and (4.2)) - // and phase for the orthogonal component of the photon field (eq. (4.3)) from the report - // "Axion-photon - // conversion in the external magnetic fields" - mpfr::mpreal CommonPhase = - Ea * 1000.0 - (axionMass * axionMass + photonMass * photonMass) / (4. * Ea * 1000.0); // in eV - mpfr::mpreal OrthogonalPhase = Ea * 1000.0 - (photonMass * photonMass) / (2. * Ea * 1000.0); // in eV - - RESTDebug << "Calculating amplitudes for one segment of trajectory where B is not zero. Segment " - "boundaries are : (" - << boundaries[i][0].X() << "," << boundaries[i][0].Y() << "," << boundaries[i][0].Z() - << ") to (" << boundaries[i][1].X() << "," << boundaries[i][1].Y() << "," - << boundaries[i][1].Z() << ")" << RESTendl; - CalculateAmplitudesInSegment(faxionAmplitude, fparallelPhotonAmplitude, forthogonalPhotonAmplitude, - averageBT, axionMass, photonMass, Ea, boundaries[i][0], boundaries[i][1], - CommonPhase, OrthogonalPhase); + 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(); - RESTDebug << RESTendl << RESTendl << RESTendl << RESTendl; - RESTDebug << "----------------------------------------------------------------" << RESTendl; - RESTDebug << " Amplitude values after calculations in one segment: " << RESTendl; - RESTDebug << " axion amplitude : "; - PrintComplex(faxionAmplitude); - RESTDebug << " parallel photon component amplitude : "; - PrintComplex(fparallelPhotonAmplitude); - RESTDebug << " orthogonal photon component amplitude : "; - PrintComplex(forthogonalPhotonAmplitude); - RESTDebug << "+--------------------------------------------------------------------------+" - << RESTendl << RESTendl; - if ((i + 1) < NofVolumes) { - cout << "Calculating amplitudes along the part of trajectory where B = 0 between the two " - "segments. Boundaries are : (" - << boundaries[i][1].X() << "," << boundaries[i][1].Y() << "," << boundaries[i][1].Z() - << ") to (" << boundaries[i + 1][0].X() << "," << boundaries[i + 1][0].Y() << "," - << boundaries[i + 1][0].Z() << ")" << RESTendl; - PropagateWithoutBField(faxionAmplitude, fparallelPhotonAmplitude, forthogonalPhotonAmplitude, - axionMass, photonMass, Ea, boundaries[i][1], boundaries[i + 1][0]); - } + fAxionEvent->SetBSquared(B * B); + fAxionEvent->SetLConversion(lenght); } - RESTDebug << RESTendl << RESTendl << RESTendl << RESTendl; - RESTDebug << "----------------------------------------------------------------" << RESTendl; - RESTDebug << " FINAL AMPLITUDES: " << RESTendl; - RESTDebug << " axion amplitude : "; - PrintComplex(faxionAmplitude); - RESTDebug << " paralell photon component amplitude : "; - PrintComplex(fparallelPhotonAmplitude); - RESTDebug << " orthogonal photon component amplitude : "; - PrintComplex(forthogonalPhotonAmplitude); - RESTDebug << " PROBABILITY for axion to photon conversion (1-|a|^2): " << 1.0 - Norm2(faxionAmplitude) - << RESTendl; - RESTDebug << "+--------------------------------------------------------------------------+" << RESTendl; - - mpfr::mpreal probabilityHighPrecision = 1.0 - Norm2(faxionAmplitude); - probability = probabilityHighPrecision.toDouble(); - fAxionEvent->SetGammaProbability(probability); - RESTDebug << "+------------------------+" << RESTendl; - RESTDebug << "Conversion probability : " << RESTendl; - RESTDebug << "fAxionEvent->GetGammaProbability() = " << fAxionEvent->GetGammaProbability() << RESTendl; - RESTDebug << "+------------------------+" << RESTendl; - - // if (fMode == "plan") - // fAxionEvent->SetPosition(MoveToPlane(position, direction, fFinalNormalPlan, fFinalPositionPlan)); - // if (fMode == "distance") fAxionEvent->SetPosition(MoveByDistance(position, direction, fDistance)); - - RESTDebug << "+------------------------+" << RESTendl; - RESTDebug << "Final position of the axion : " << RESTendl; - RESTDebug << "(" << fAxionEvent->GetPositionX() << "," << fAxionEvent->GetPositionY() << "," - << fAxionEvent->GetPositionZ() << ")" << RESTendl; - RESTDebug << "+------------------------+" << RESTendl; if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fAxionEvent->PrintEvent(); - boundaries.clear(); + /// Missing to propagate the axion to the end of magnet bore? return fAxionEvent; } - -/////////////////////////////////////////////// -/// \brief Function reading input parameters from the RML TRestAxionFieldPropagationProcess metadata section -/// -void TRestAxionFieldPropagationProcess::InitFromConfigFile() { - this->Initialize(); - - fMode = GetParameter("mode"); - fFinalNormalPlan = Get3DVectorParameterWithUnits("finalNPlan"); - fFinalPositionPlan = Get3DVectorParameterWithUnits("finalPositionPlan"); - fDistance = GetDblParameterWithUnits("distance"); - - PrintMetadata(); -} diff --git a/src/TRestAxionGeneratorProcess.cxx b/src/TRestAxionGeneratorProcess.cxx index c53de957..05216662 100644 --- a/src/TRestAxionGeneratorProcess.cxx +++ b/src/TRestAxionGeneratorProcess.cxx @@ -21,16 +21,41 @@ *************************************************************************/ ////////////////////////////////////////////////////////////////////////// -/// TRestAxionGeneratorProcess TOBE documented /// -/// The axion is generated with intensity proportional to g_ag = 1.0 x g10 +/// This generator will produce solar axion events with the Sun at position +/// (0,0,-AU) and detector at position (0,0,0). +/// +/// The axion generation properties, such as coupling type and intensity, are +/// defined by TRestAxionSolarFlux. That class defines the method +/// TRestAxionSolarFlux::GetRandomEnergyAndRadius that is exploited by this +/// process to define initial conditions for each axion, such as angular +/// direction or energy. The flux intensity and its dependency with the solar +/// radius and spectral energy will be fully coded inside TRestAxionSolarFlux, +/// while this method will be able to play with the Sun-generator position, +/// and the size of the target detector (which is placed at (0,0,0)) to +/// define an initial axion position, direction and energy. +/// +/// The following metadata members can be tuned in this process: +/// - **axionMass**: It defines the axion mass in eV required later on by axion-photon +/// conversion processes. Its default value is 0. +/// - **targetRadius**: All generated axion events will end up in a circular region +/// placed at the XY plane. This parameter defines the radius of the circular region. +/// The default value is 80cm. It will be used to generate an additional deviation +/// to the position. +/// - **generatorType**: It defines for the moment two different types of generator. +/// - `solarFlux`: It places the particle at the solar disk, 1-AU distance, with +/// a spatial and energy distribution given by the TRestAxionSolarFlux description. +/// - `flat`: It just gives a parallel flux with the extension of fTargetRadius. +/// The photons energy is fixed for the moment to 1keV, and the Z-position far enough +/// to be outside the IAXO helioscope. +/// ///-------------------------------------------------------------------------- /// /// RESTsoft - Software for Rare Event Searches with TPCs /// /// History of developments: /// -/// 2019-March: First implementation of shared memory buffer to rawsignal conversion. +/// 2022-February: New generator implementation using TRestAxionSolarFlux /// Javier Galan /// /// \class TRestAxionGeneratorProcess @@ -102,8 +127,7 @@ void TRestAxionGeneratorProcess::Initialize() { fOutputAxionEvent = new TRestAxionEvent(); fIsExternal = true; - - fRandom = new TRandom3(0); + fSingleThreadOnly = true; } /////////////////////////////////////////////// @@ -113,184 +137,98 @@ void TRestAxionGeneratorProcess::Initialize() { void TRestAxionGeneratorProcess::InitProcess() { RESTDebug << "Entering ... TRestAxionGeneratorProcess::InitProcess" << RESTendl; - fAxionSpectrum = (TRestAxionSpectrum*)this->GetMetadata("TRestAxionSpectrum"); + fAxionFlux = GetMetadata(); - if (!fAxionSpectrum) { - RESTError << "TRestAxionGeneratorProcess. Axion model was not defined!" << RESTendl; - exit(0); + if (fGeneratorType == "solarFlux" && fAxionFlux == nullptr) { + if (!this->GetError()) this->SetError("The solar flux definition was not found."); } -} - -/////////////////////////////////////////////// -/// \brief This method gets a random energy relying on the solar axion model defined in TRestAxionSolarModel -/// -Double_t TRestAxionGeneratorProcess::GenerateEnergy() { - RESTDebug << "Entering TRestAxionGeneratorProcess::GenerateEnergy() ..." << RESTendl; - Double_t solarFlux = fAxionSpectrum->GetSolarAxionFlux(fEnergyRange.X(), fEnergyRange.Y(), fEnergyStep); - - Double_t random = solarFlux * fRandom->Uniform(0, 1.0); - - Double_t sum = 0; - for (double en = fEnergyRange.X(); en < fEnergyRange.Y(); en += fEnergyStep) { - sum += fAxionSpectrum->GetDifferentialSolarAxionFlux(en) * fEnergyStep; - if (random < sum) { - RESTDebug << "TRestAxionGeneratorProcess::GenerateEnergy. Energy = " << en << RESTendl; - return en + fRandom->Uniform(0, fEnergyStep); - } + if (!fRandom) { + delete fRandom; + fRandom = nullptr; } - return fEnergyRange.Y(); + fRandom = new TRandom3(fSeed); + if (fSeed == 0) fSeed = fRandom->GetSeed(); } -TVector3 TRestAxionGeneratorProcess::GenerateDirection() { - TVector3 direction; - - // For the moment "circleWall" just generates in the XY-plane - if (fAngularDistribution == "flux") { - return fAngularDirection; - } - - RESTWarning << "Angular distribution : " << fAngularDistribution << " is not defined!" << RESTendl; - - return fAngularDirection; -} - -TVector3 TRestAxionGeneratorProcess::GeneratePosition() { - TVector3 position; - Double_t r, x, y, z; - - if (fSpatialDistribution == "circleWall") { - do { - x = fRandom->Uniform(-fSpatialRadius, fSpatialRadius); - y = fRandom->Uniform(-fSpatialRadius, fSpatialRadius); - - r = x * x + y * y; - } while (r > fSpatialRadius * fSpatialRadius); - - position = TVector3(x, y, z); //+ fSpatialOrigin; - - if (fMode == "XYZRotation") { - position.RotateX(fRotation.X()); - position.RotateY(fRotation.Y()); - position.RotateZ(fRotation.Z()); - position = position + fSpatialOrigin; - } - - if (fMode == "nPlanRotation") { - fNormalPlan = fNormalPlan.Unit(); - TVector3 rotVec(0, -fNormalPlan.X(), fNormalPlan.Y()); - Double_t angle = fNormalPlan.Angle(TVector3(1, 0, 0)); - position.Rotate(angle, rotVec); - position = position + fSpatialOrigin; - } - - return position; - } - - if (fSpatialDistribution == "circleWallXY") { - do { - x = fRandom->Uniform(-fSpatialRadius, fSpatialRadius); - y = fRandom->Uniform(-fSpatialRadius, fSpatialRadius); - - r = x * x + y * y; - } while (r > fSpatialRadius * fSpatialRadius); - - position = TVector3(x, y, 0) + fSpatialOrigin; - - return position; - } - - if (fSpatialDistribution == "circleWallXZ") { - do { - x = fRandom->Uniform(-fSpatialRadius, fSpatialRadius); - z = fRandom->Uniform(-fSpatialRadius, fSpatialRadius); - - r = x * x + z * z; - } while (r > fSpatialRadius * fSpatialRadius); - - position = TVector3(x, 0, z) + fSpatialOrigin; - - return position; - } - - if (fSpatialDistribution == "circleWallYZ") { - do { - y = fRandom->Uniform(-fSpatialRadius, fSpatialRadius); - z = fRandom->Uniform(-fSpatialRadius, fSpatialRadius); +/////////////////////////////////////////////// +/// \brief The main processing event function +/// +TRestEvent* TRestAxionGeneratorProcess::ProcessEvent(TRestEvent* evInput) { + RESTDebug << "TRestAxionGeneratorProcess::ProcessEvent : " << fCounter << RESTendl; + fOutputAxionEvent->SetID(fCounter); + fCounter++; - r = y * y + z * z; - } while (r > fSpatialRadius * fSpatialRadius); + TVector3 axionPosition = TVector3(0, 0, -REST_Physics::AU); + TVector3 axionDirection = TVector3(0, 0, 1); + Double_t energy = 1; - position = TVector3(0, y, z) + fSpatialOrigin; + // Random unit position + // We avoid the use of expensive trigonometric functions + Double_t x, y, r; + do { + x = 2 * (fRandom->Rndm() - 0.5); + y = 2 * (fRandom->Rndm() - 0.5); + r = x * x + y * y; + } while (r > 1 || r == 0); - return position; - } + r = TMath::Sqrt(r); - if (fSpatialDistribution == "sphereIn") { - fRandom->Sphere(x, y, z, fSpatialRadius); - position = TVector3(x, y, z) + fSpatialOrigin; + if (fGeneratorType == "solarFlux") { + std::pair p = fAxionFlux->GetRandomEnergyAndRadius(); + energy = p.first; + Double_t radius = p.second; - fAngularDirection = -TVector3(x, y, z); - fAngularDirection = fAngularDirection.Unit(); + axionPosition = TVector3(REST_Physics::solarRadius * radius * x, + REST_Physics::solarRadius * radius * y, -REST_Physics::AU); - return position; + axionDirection = -axionPosition.Unit(); } - if (fSpatialDistribution == "sphereOut") { - fRandom->Sphere(x, y, z, fSpatialRadius); - position = TVector3(x, y, z) + fSpatialOrigin; - - fAngularDirection = TVector3(x, y, z); - fAngularDirection = fAngularDirection.Unit(); - - return position; + if (fGeneratorType == "flat" || fGeneratorType == "plain") { + if (fMaxEnergy > 0) + energy = fRandom->Rndm() * (fMaxEnergy - fMinEnergy) + fMinEnergy; + else + energy = (1. - fMinEnergy) * fRandom->Rndm() + fMinEnergy; } - RESTWarning << "Spatial distribution : " << fSpatialDistribution << " is not defined!" << RESTendl; - - return position; -} - -/////////////////////////////////////////////// -/// \brief The main processing event function -/// -TRestEvent* TRestAxionGeneratorProcess::ProcessEvent(TRestEvent* evInput) { - RESTDebug << "TRestAxionGeneratorProcess::ProcessEvent : " << fCounter << RESTendl; - fOutputAxionEvent->SetID(fCounter); - fCounter++; - - fOutputAxionEvent->SetEnergy(GenerateEnergy()); - fOutputAxionEvent->SetPosition(GeneratePosition()); - fOutputAxionEvent->SetDirection(GenerateDirection()); + /// The axion position must be displaced by the target size. + /// We always do this. It is independent of generator + /// The target is virtually placed at the (0,0,0). + /// In my opinion the target should be either the optics, or the magnet end bore. + /// Then one should place the optics or the magnet end bore at the (0,0,0). + /// + do { + x = 2 * (fRandom->Rndm() - 0.5); + y = 2 * (fRandom->Rndm() - 0.5); + r = x * x + y * y; + } while (r > 1 || r == 0); + + r = TMath::Sqrt(r); + + axionPosition = axionPosition + TVector3(fTargetRadius * x, fTargetRadius * y, 0); + /// /// + + fOutputAxionEvent->SetEnergy(energy); + fOutputAxionEvent->SetPosition(axionPosition); + fOutputAxionEvent->SetDirection(axionDirection); fOutputAxionEvent->SetMass(fAxionMass); - // cout << "anaTree : " << fAnalysisTree << endl; - // fAnalysisTree->SetObservableValue( this, "energy", fOutputAxionEvent->GetEnergy() ); - // fAnalysisTree->PrintObservables(); - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fOutputAxionEvent->PrintEvent(); return fOutputAxionEvent; } -/////////////////////////////////////////////// -/// \brief Function reading input parameters from the RML TRestAxionGeneratorProcess metadata section -/// -void TRestAxionGeneratorProcess::InitFromConfigFile() { - // These 2-params should be moved to TRestAxionSolarModel - fEnergyRange = Get2DVectorParameterWithUnits("energyRange"); - fEnergyStep = GetDblParameterWithUnits("energyStep"); - - fAngularDistribution = GetParameter("angularDistribution", "flux"); - fAngularDirection = Get3DVectorParameterWithUnits("angularDirection"); +void TRestAxionGeneratorProcess::PrintMetadata() { + TRestMetadata::PrintMetadata(); - fSpatialDistribution = GetParameter("spatialDistribution", "circleWall"); - fSpatialRadius = GetDblParameterWithUnits("spatialRadius"); - fSpatialOrigin = Get3DVectorParameterWithUnits("spatialOrigin"); + RESTMetadata << "Generator type: " << fGeneratorType << RESTendl; + RESTMetadata << "Axion mass: " << fAxionMass * units("eV") << " eV" << RESTendl; + RESTMetadata << "Target radius: " << fTargetRadius * units("cm") << " cm" << RESTendl; + RESTMetadata << "Random seed: " << (UInt_t)fSeed << RESTendl; + RESTMetadata << "Energy range: (" << fMinEnergy << ", " << fMaxEnergy << ") keV" << RESTendl; - fRotation = StringTo3DVector(GetParameter("rotation")); - fNormalPlan = Get3DVectorParameterWithUnits("normalWall"); - fMode = GetParameter("mode"); + RESTMetadata << "+++++++++++++++++++++++++++++++++++++++++++++++++" << RESTendl; } diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx index 50d34c8f..4e68edcf 100644 --- a/src/TRestAxionMagneticField.cxx +++ b/src/TRestAxionMagneticField.cxx @@ -1154,13 +1154,11 @@ std::vector TRestAxionMagneticField::GetTransversalComponentAlongPath( /// Double_t TRestAxionMagneticField::GetTransversalFieldAverage(TVector3 from, TVector3 to, Double_t dl, Int_t Nmax) { - Double_t length = (to - from).Mag(); - Double_t Bavg = 0.; std::vector Bt = GetTransversalComponentAlongPath(from, to, dl, Nmax); for (auto& b : Bt) Bavg += b; - if (length > 0) return Bavg / length; + if (Bt.size() > 0) return Bavg / Bt.size(); RESTError << "TRestAxionMagneticField::GetTransversalFieldAverage. Lenght is zero!" << RESTendl; return 0.; diff --git a/src/TRestAxionOptics.cxx b/src/TRestAxionOptics.cxx index 401fd6f6..40382aed 100644 --- a/src/TRestAxionOptics.cxx +++ b/src/TRestAxionOptics.cxx @@ -344,6 +344,30 @@ void TRestAxionOptics::ResetPositions() { fCurrentMirror = -1; } +/////////////////////////////////////////////// +/// \brief It returns the last valid particle position known in the particle tracking. +/// +TVector3 TRestAxionOptics::GetLastGoodPosition() { + if (GetExitDirection().Mag() > 0) + return GetExitPosition(); + else if (GetMiddleDirection().Mag() > 0) + return GetMiddlePosition(); + + return GetEntrancePosition(); +} + +/////////////////////////////////////////////// +/// \brief It returns the last valid particle direction known in the particle tracking. +/// +TVector3 TRestAxionOptics::GetLastGoodDirection() { + if (GetExitDirection().Mag() > 0) + return GetExitDirection(); + else if (GetMiddleDirection().Mag() > 0) + return GetMiddleDirection(); + + return GetEntranceDirection(); +} + /////////////////////////////////////////////// /// \brief Propagating photon /// @@ -416,6 +440,8 @@ void TRestAxionOptics::PrintMetadata() { RESTMetadata << "---------" << RESTendl; RESTMetadata << "Entrance position in Z : " << GetEntrancePositionZ() << " mm" << RESTendl; RESTMetadata << "Exit position in Z : " << GetExitPositionZ() << " mm" << RESTendl; + RESTMetadata << "Low radial limit : " << GetRadialLimits().first << " mm" << RESTendl; + RESTMetadata << "High radial limit : " << GetRadialLimits().second << " mm" << RESTendl; RESTMetadata << "---------" << RESTendl; RESTMetadata << " " << RESTendl; RESTMetadata << " Use \"this->PrintMasks()\" to get masks info" << RESTendl; diff --git a/src/TRestAxionOpticsProcess.cxx b/src/TRestAxionOpticsProcess.cxx new file mode 100644 index 00000000..d72e9ab8 --- /dev/null +++ b/src/TRestAxionOpticsProcess.cxx @@ -0,0 +1,161 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TRestAxionOpticsProcess TOBE documented +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2019-March: First implementation of a dummy optics response +/// Javier Galan +/// +/// \class TRestAxionOpticsProcess +/// \author +/// +///
+/// +#include "TRestAxionOpticsProcess.h" +using namespace std; + +ClassImp(TRestAxionOpticsProcess); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestAxionOpticsProcess::TRestAxionOpticsProcess() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// +TRestAxionOpticsProcess::TRestAxionOpticsProcess(char* cfgFileName) { + Initialize(); + + LoadConfig(cfgFileName); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestAxionOpticsProcess::~TRestAxionOpticsProcess() { delete fAxionEvent; } + +/////////////////////////////////////////////// +/// \brief Function to load the default config in absence of RML input +/// +void TRestAxionOpticsProcess::LoadDefaultConfig() { + SetName(this->ClassName()); + SetTitle("Default config"); +} + +/////////////////////////////////////////////// +/// \brief Function to load the configuration from an external configuration file. +/// +/// If no configuration path is defined in TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// correspondig TRestGeant4AnalysisProcess section inside the RML. +/// +void TRestAxionOpticsProcess::LoadConfig(std::string cfgFilename, std::string name) { + if (LoadConfigFromFile(cfgFilename, name)) LoadDefaultConfig(); +} + +/////////////////////////////////////////////// +/// \brief Function to initialize input/output event members and define the section name +/// +void TRestAxionOpticsProcess::Initialize() { + SetSectionName(this->ClassName()); + SetLibraryVersion(LIBRARY_VERSION); + + fAxionEvent = new TRestAxionEvent(); +} + +/////////////////////////////////////////////// +/// \brief Process initialization. Data members that require initialization just before start processing +/// should be initialized here. +/// +void TRestAxionOpticsProcess::InitProcess() { + RESTDebug << "Entering ... TRestAxionGeneratorProcess::InitProcess" << RESTendl; + + fOptics = GetMetadata(); + + if (fOptics) { + fOptics->PrintMetadata(); + } else { + RESTError << "TRestAxionOptics::InitProcess. No sucess instantiating optics." << RESTendl; + } +} + +/////////////////////////////////////////////// +/// \brief The main processing event function +/// +TRestEvent* TRestAxionOpticsProcess::ProcessEvent(TRestEvent* evInput) { + fAxionEvent = (TRestAxionEvent*)evInput; + + TVector3 inPos = fAxionEvent->GetPosition(); + TVector3 inDir = fAxionEvent->GetDirection(); + Double_t energy = fAxionEvent->GetEnergy(); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fAxionEvent->PrintEvent(); + Double_t efficiency = fOptics->PropagatePhoton(inPos, inDir, energy); + + SetObservableValue("efficiency", efficiency); + + // We register the event even if it is not properly reflected, i.e. efficiency = 0 + fAxionEvent->SetPosition(fOptics->GetLastGoodPosition()); + fAxionEvent->SetDirection(fOptics->GetLastGoodDirection()); + + /// TODO set efficiency inside the event + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { + fAxionEvent->PrintEvent(); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) GetChar(); + } + + 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 new file mode 100644 index 00000000..ef901b19 --- /dev/null +++ b/src/TRestAxionTransportProcess.cxx @@ -0,0 +1,149 @@ +/************************************************************************* + * 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. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TRestAxionTransportProcess simply translates the axion to a given +/// z-position without modifying the direction. This process will only +/// allow the movement on the sense of direction, particle cannot move +/// backwards. +/// +/// It receives a single parameter that defines the Z-position where +/// the particle should be moved `zPosition`. +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2022-June: Basic implementation of a photon transport process +/// Javier Galan +/// +/// \class TRestAxionTransportProcess +/// \author +/// +///
+/// +#include "TRestAxionTransportProcess.h" +using namespace std; + +ClassImp(TRestAxionTransportProcess); + +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestAxionTransportProcess::TRestAxionTransportProcess() { Initialize(); } + +/////////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// +TRestAxionTransportProcess::TRestAxionTransportProcess(char* cfgFileName) { + Initialize(); + + LoadConfig(cfgFileName); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestAxionTransportProcess::~TRestAxionTransportProcess() { delete fAxionEvent; } + +/////////////////////////////////////////////// +/// \brief Function to load the default config in absence of RML input +/// +void TRestAxionTransportProcess::LoadDefaultConfig() { + SetName(this->ClassName()); + SetTitle("Default config"); +} + +/////////////////////////////////////////////// +/// \brief Function to load the configuration from an external configuration file. +/// +/// If no configuration path is defined in TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// correspondig TRestGeant4AnalysisProcess section inside the RML. +/// +void TRestAxionTransportProcess::LoadConfig(std::string cfgFilename, std::string name) { + if (LoadConfigFromFile(cfgFilename, name)) LoadDefaultConfig(); +} + +/////////////////////////////////////////////// +/// \brief Function to initialize input/output event members and define the section name +/// +void TRestAxionTransportProcess::Initialize() { + SetSectionName(this->ClassName()); + SetLibraryVersion(LIBRARY_VERSION); + + fAxionEvent = new TRestAxionEvent(); +} + +/////////////////////////////////////////////// +/// \brief Process initialization. Data members that require initialization just before start processing +/// should be initialized here. +/// +void TRestAxionTransportProcess::InitProcess() { + RESTDebug << "Entering ... TRestAxionGeneratorProcess::InitProcess" << RESTendl; +} + +/////////////////////////////////////////////// +/// \brief The main processing event function +/// +TRestEvent* TRestAxionTransportProcess::ProcessEvent(TRestEvent* evInput) { + fAxionEvent = (TRestAxionEvent*)evInput; + + TVector3 inPos = fAxionEvent->GetPosition(); + TVector3 inDir = fAxionEvent->GetDirection(); + + if ((inPos.Z() - fZPosition) * inDir.Z() >= 0) { + RESTWarning + << "TRestAxionTransportProcess::ProcessEvent. Not appropiate particle direction to reach Z=" + << fZPosition << RESTendl; + RESTWarning << "Axion position. Z:" << inPos.Z() << " direction Z: " << inDir.Z() << RESTendl; + fAxionEvent->PrintEvent(); + + return fAxionEvent; + } + + TVector3 newPosition = + REST_Physics::MoveToPlane(inPos, inDir, TVector3(0, 0, 1), TVector3(0, 0, fZPosition)); + + fAxionEvent->SetPosition(newPosition); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) { + fAxionEvent->PrintEvent(); + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) GetChar(); + } + + return fAxionEvent; +}