rest-for-physics · jgalan · Feb 9, 2022 · Feb 9, 2022 · Feb 9, 2022 · Feb 9, 2022
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -61,6 +61,11 @@ endif (${REST_SOLAXFLUX} MATCHES "ON")
 
 set( external_libs ${external_libs} -lmpfr )
 
+file(GLOB_RECURSE addon_src
+    "mcpl.c" )
+
+set(addon_inc ${CMAKE_CURRENT_SOURCE_DIR}/external/mcpl)
+
 COMPILELIB("")
 
 INSTALL(DIRECTORY ./data/

diff --git 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](https://doi.org/10.1016/j.cpc.2017.04.012).
diff --git a/external/mcpl/mcpl.h b/external/mcpl/mcpl.h
diff --git a/inc/TRestAxionGeneratorProcess.h b/inc/TRestAxionGeneratorProcess.h
@@ -40,7 +40,7 @@ class TRestAxionGeneratorProcess : public TRestEventProcess {
     Int_t fCounter = 0;  //!
 
     /// A pointer to the axion model stored in TRestRun
-    TRestAxionSpectrum *fAxionSpectrum; //!
+    TRestAxionSpectrum* fAxionSpectrum;  //!
 
     /// Random number generator
     TRandom3* fRandom;  //!

diff --git a/inc/TRestAxionLikelihood.h b/inc/TRestAxionLikelihood.h
@@ -25,10 +25,10 @@
 
 #include <TRestMetadata.h>
 
-#include "TRestAxionSpectrum.h"
 #include "TRestAxionBufferGas.h"
 #include "TRestAxionPhotonConversion.h"
 #include "TRestAxionSolarModel.h"
+#include "TRestAxionSpectrum.h"
 
 #include "TRandom3.h"
 
@@ -63,9 +63,9 @@ class TRestAxionLikelihood : public TRestMetadata {
 
     Double_t fLastStepDensity = 0.;  //->
 
-    TRestAxionPhotonConversion *fPhotonConversion; //!
-    TRestAxionBufferGas *fBufferGas; //!
-    TRestAxionSpectrum *fAxionSpectrum; //!
+    TRestAxionPhotonConversion* fPhotonConversion;  //!
+    TRestAxionBufferGas* fBufferGas;                //!
+    TRestAxionSpectrum* fAxionSpectrum;             //!
 
     /// Random number generator
     TRandom3* fRandom;  //!

diff --git a/inc/TRestAxionMCPLOptics.h b/inc/TRestAxionMCPLOptics.h
@@ -38,6 +38,8 @@ class TRestAxionMCPLOptics : public TRestAxionOptics {
     std::string fOutputMCPLFilename;
 
    public:
+    Int_t LoadMCPLFiles();
+
     void PrintMetadata();
 
     void InitFromConfigFile();

diff --git a/inc/mpreal.h b/inc/mpreal.h
@@ -125,9 +125,10 @@
 #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)
+#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

diff --git a/pipeline/optics/setups.rml b/pipeline/optics/setups.rml
@@ -1,5 +1,5 @@
 <xml>
-<TRestAxionMCPLOptics name="mcpl" >
+<TRestAxionMCPLOptics name="mcpl" verboseLevel="info" >
 	<parameter name="center" value="(0,0,200)mm" />
 	<parameter name="axis" value="(0,0.02,0.98)" />
 	<parameter name="length" value="22cm" />
@@ -8,9 +8,9 @@
 	<parameter name="shellMinRadii" value="5,10,15,20,25" />
 	<parameter name="shellMaxRadii" value="9.9,14.9,19.9,24.9,29.9" />
 
-	<parameter name="inputMCPLFilename" value="dummyInput.mcpl" />
-	<parameter name="outputMCPLFilename" value="dummyOutput.mcpl" />
 
+	<parameter name="inputMCPLFilename" value="https://sultan.unizar.es/mcpl/mcpl_detector_plane.mcpl" />
+	<parameter name="outputMCPLFilename" value="https://sultan.unizar.es/mcpl/mcpl_post_optics.mcpl" />
 </TRestAxionMCPLOptics>
 
 <TRestAxionOptics name="mcpl2" center="(0,0,300)cm" axis="(0,0,3)" />

diff --git a/src/TRestAxionGeneratorProcess.cxx b/src/TRestAxionGeneratorProcess.cxx
@@ -126,8 +126,7 @@ void TRestAxionGeneratorProcess::InitProcess() {
 ///
 Double_t TRestAxionGeneratorProcess::GenerateEnergy() {
     debug << "Entering TRestAxionGeneratorProcess::GenerateEnergy() ..." << endl;
-    Double_t solarFlux =
-        fAxionSpectrum->GetSolarAxionFlux(fEnergyRange.X(), fEnergyRange.Y(), fEnergyStep);
+    Double_t solarFlux = fAxionSpectrum->GetSolarAxionFlux(fEnergyRange.X(), fEnergyRange.Y(), fEnergyStep);
 
     Double_t random = solarFlux * fRandom->Uniform(0, 1.0);
 

diff --git a/src/TRestAxionMCPLOptics.cxx b/src/TRestAxionMCPLOptics.cxx
@@ -63,6 +63,7 @@
 ///
 
 #include "TRestAxionMCPLOptics.h"
+#include "mcpl.h"
 
 using namespace std;
 
@@ -114,6 +115,54 @@ void TRestAxionMCPLOptics::Initialize() {
 
     SetSectionName(this->ClassName());
     SetLibraryVersion(LIBRARY_VERSION);
+
+    if (fInputMCPLFilename != "" && fOutputMCPLFilename != "") LoadMCPLFiles();
+}
+
+///////////////////////////////////////////////
+/// \brief Method to load inside the class the MCPL files
+///
+Int_t TRestAxionMCPLOptics::LoadMCPLFiles() {
+    string fname = fInputMCPLFilename;
+    if (fname.find("https") == 0) fname = TRestTools::DownloadRemoteFile(fname);
+    mcpl_file_t f = mcpl_open_file(fname.c_str());
+
+    string fname2 = fOutputMCPLFilename;
+    if (fname2.find("https") == 0) fname2 = TRestTools::DownloadRemoteFile(fname2);
+    mcpl_file_t g = mcpl_open_file(fname2.c_str());
+
+    const mcpl_particle_t *p, *q;
+    while ((p = mcpl_read(f))) {
+        q = mcpl_read(g);
+
+        cout << "Input MCPL file (detector_plane)" << endl;
+        cout << "--------------------------------" << endl;
+        cout << "X: " << p->position[0] << " Y: " << p->position[1] << " Z: " << p->position[2] << endl;
+        cout << "pX: " << p->direction[0] << " pY: " << p->direction[1] << " pZ: " << p->direction[2] << endl;
+        cout << "polX: " << p->polarisation[0] << " polY: " << p->polarisation[1]
+             << " polZ: " << p->polarisation[2] << endl;
+        cout << "Ekin = " << p->ekin << endl;
+        cout << "Time = " << p->time << endl;
+        cout << "Weight = " << p->weight << endl;
+        cout << "User flags = " << p->userflags << endl;
+        cout << endl;
+        cout << "Output MCPL file (post_optics)" << endl;
+        cout << "--------------------------------" << endl;
+        cout << "X: " << q->position[0] << " Y: " << q->position[1] << " Z: " << q->position[2] << endl;
+        cout << "pX: " << q->direction[0] << " pY: " << q->direction[1] << " pZ: " << q->direction[2] << endl;
+        cout << "polX: " << q->polarisation[0] << " polY: " << q->polarisation[1]
+             << " polZ: " << q->polarisation[2] << endl;
+        cout << "Ekin = " << q->ekin << endl;
+        cout << "Time = " << q->time << endl;
+        cout << "Weight = " << q->weight << endl;
+        cout << "User flags = " << q->userflags << endl;
+        cout << endl;
+
+        getchar();
+    }
+    mcpl_close_file(f);
+
+    return 0;
 }
 
 ///////////////////////////////////////////////

diff --git a/src/TRestAxionMagneticField.cxx b/src/TRestAxionMagneticField.cxx
@@ -1080,9 +1080,9 @@ TVector3 TRestAxionMagneticField::GetFieldAverageTransverseVector(TVector3 from,
 
     if ((length > 0) && (numberofpoints > 0)) {
         Bavg = Bavg * (1.0 / numberofpoints);  // calculates the average magnetic field vector
-        BTavg = Bavg -
-                (Bavg * direction) *
-                    direction;  // calculates the transverse component of the average magnetic field vector
+        BTavg =
+            Bavg - (Bavg * direction) *
+                       direction;  // calculates the transverse component of the average magnetic field vector
         debug << "B average vector = (" << Bavg.x() << ", " << Bavg.y() << ", " << Bavg.z() << ")" << endl;
         debug << "Transverse B average vector = (" << BTavg.x() << ", " << BTavg.y() << ", " << BTavg.z()
               << ")" << endl;

diff --git a/src/TRestAxionSolarModel.cxx b/src/TRestAxionSolarModel.cxx
@@ -20,7 +20,7 @@
  * For the list of contributors see $REST_PATH/CREDITS.                  *
  *************************************************************************/
 
-/***************** DOXYGEN DOCUMENTATION ********************************
+//////////////////////////////////////////////////////////////////////////
 /// TRestAxionSolarModel is a class used to calculate all axion interaction
 /// rates in the Sun based on a solar model and opacity code. The main
 /// purpose of the class is to return the two types of solar axion spectra,
@@ -49,11 +49,11 @@
 ///             Sebastian Hoof
 ///
 /// \class      TRestAxionSolarModel
-/// \author     Javier Galan, Sebastian Hoof
+/// \author     Javier Galan
+/// \author	    Sebastian Hoof
 ///
 /// <hr>
 ///
- *************************************************************************/
 
 #include "TRestAxionSolarModel.h"