Skip to content

position dependent photon field for EM interactions (optimised)#535

Open
GDMarco wants to merge 120 commits intoCRPropa:masterfrom
GDMarco:PositionDependentPhotonFieldOpt
Open

position dependent photon field for EM interactions (optimised)#535
GDMarco wants to merge 120 commits intoCRPropa:masterfrom
GDMarco:PositionDependentPhotonFieldOpt

Conversation

@GDMarco
Copy link

@GDMarco GDMarco commented Jan 26, 2026

General description

This is a fully working version of the CRPropa code with some implementation to allow spatially dependent photon field to be initialised and employed for electromagnetic interactions. A particular example is the introduction of the galactic interstellar field (the density and interaction rate tables are at this repo: ISRF tables.)
In principle, this implementation can be used for any spatial dependent photon field properly build by the user.
To allow the introduction of such spatially variable photon fields in CRPropa, the classe and modules that have been modified arePhotonBackground.*, EMPairProduction.*, EMDoublePairProduction.*, EMInverseComptonScattering.* and EMTripletPairProduction.*. Moreover, the class InteractionRates.* has been implemented to handle the process rates.
An example on how to use the code is contained in the notebooks.

Technical details

The .cc files of PhotonBackground.*,EMPairProduction.*, EMDoublePairProduction.*, EMInverseComptonScattering.* and EMTripletPairProduction.* employ some tools of the #include <filesystem> library, only available since C++17. The proper library is specified in the CMakeLists.txt files, as well as in the header of the .cc files.

The code uses the nanoflann functionalities to look for the closest photon field point in your custom grid and evaluate the density and/or the rate correctly. In this new version, nanoflann is implemented as a header-only library among the ones in the /libs/ directory.

GDMarco added 30 commits March 5, 2024 16:40
//InteractionRatesIsotropic::InteractionRatesIsotropic();
//InteractionRatesIsotropic::InteractionRatesPositionDependent();
Changing nomenclature (isSpatialDependent -> isPositionDependent).
Changing nomenclature
Change nomenclature.
Change nomenclature.
style & constructor of InteractionRates subclasses
@GDMarco GDMarco marked this pull request as ready for review February 12, 2026 14:24
Copy link
Contributor

@JanNiklasB JanNiklasB left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi,

thanks for your contribution!
I have done a first review of your pull request and noticed the following things.

  • Please try to not include lines in your pull request that are only altered in white spaces without having any purpose, you often only change tabs to spaces or add new empty lines, this makes it difficult to spot actual changes and leads to difficult to read code.
  • Reduce the amount of includes you use in your files, any not needed include or already in the header existing includes should not be in the .cpp file again. This is only a style thing since otherwise the code seems a little cluttered.
  • For the 4 changed module files (EM*.h and EM*.cpp) I only commented in EMDoublePairProduction.* but my comments hold true for all module files:
    -> The old behavior of function and classes regarding function calls can not be changed, it needs to be backwards compatible ( We need to at least notify the community first when we remove or change crucial features)
    -> Everything needs to be as modular as possible, so if you use a class like the interaction rates in a function you need to only require the base class so the user can create their own versions and hand them over. Those base classes then need everything to distinguish them inside those functions
  • You require the function splitFilename in many of your classes but this function does not depend on any class member, so it would be better to add this function to Common.h as a standalone
  • Every function which is not obvious to use or where it is not obvious what it does needs to have a description
  • You need to add your changes to the change log CHANGELOG.md

Many comments would repeat themself so I usually mention something only 1 or 2 times.

With regards,
Jan-Niklas

#include "crpropa/Random.h"
#include "crpropa/PhotonBackground.h"
#include "crpropa/InteractionRates.h"
#include "crpropa/Geometry.h"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PhotonBackground.h, InteractionRates.h and Geometry.h are already included in the corresponding EMDoublePairProduction.h header. Please remove them here for style reasons.

@param thinning weighted sampling of secondaries (0: all particles are tracked; 1: maximum thinning)
@param limit step size limit as fraction of mean free path
@param surface suface to enclose the grid nodes to be loaded
@param interactionRates(hidden) object to store and access to the interaction rates of the process
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The interactionRates should not appear as parameter since the user is not able to actually specify them that way. Instead mention them in the description of the constructor like for example:

/** Constructor
* You can set interactionRates via initRate(std:.string filename)
* ....
* /

* @param surface closed surface to confine the grid to be uploaded
*/
void setSurface(ref_ptr<Surface> surface);
bool hasSurface() const;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It might be better to return the Surface pointer instead, that way the user can do the check and still use the pointer afterwards.

std::string getInteractionTag() const;

void initRate(std::string filename, InteractionRatesHomogeneous* intRatesHom);
void initRatePositionDependentPhotonField(std::string filepath, InteractionRatesPositionDependent* intRatesPosDep);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A description of initRate and initRatePositionDependentPhotonField is necessary, it is not clear for the average user how to use those functions.

Furthermore, we might need to add a initRate(std::string) overload so the old call behavior is preserverd.

if (!this->photonField->hasPositionDependence()) {

this->interactionRates = new InteractionRatesHomogeneous("interactionRatesHomogeneous", false);
InteractionRatesHomogeneous* intRatesHom = static_cast<InteractionRatesHomogeneous*>(this->interactionRates.get());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With the requested change of initRate in the header file this pointer and the one in line 51 are not required. Furthermore initRate would be used in both cases so the if condition would is then only required inside the initRate function.

initRedshiftScaling();
this->isPositionDependent = isPositionDependent;

if (this->isPositionDependent) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

instead of throwing an error when the user sets isPositionDependent=true remove the argument from the constructor and just set it manually to false

} else if (!this->isPositionDependent) {

throw std::runtime_error("Photon Field " + fieldName + " is not position dependent! It is not the correct class. \n");

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Like above, if you do not want to have a RedshiftDependence then don't give the user the option to set it

*/
class TabularSpatialPhotonField: public PhotonField {
public:
TabularSpatialPhotonField(const std::string fieldName, const bool isRedshiftDependent = true, const bool isPositionDependent = true, ref_ptr<Surface> surface = nullptr);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the field is not dependent on Redshift or Position, then it would be better to manually set them to false instead of giving the user the option to set them here. Also, since the options are set to true by default and the constructor throws an error when that happens, the user probably is forced to set them to false in every case.

if(APPLE)
message(STATUS "Building on macOS, using C++11")
set(CMAKE_CXX_STANDARD 11)
# set(CMAKE_CXX_STANDARD_REQUIRED ON)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there something on macOS that prevents us to set this CXX standard to required? If not it would be good to uncomment this I think.


endif(ENABLE_PYTHON AND Python_FOUND)

target_include_directories(crpropa
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add nanoflann to crpropa over CRPROPA_EXTRA_INCLUDE and CRPROPA_EXTRA_SOURCES (see other libraries).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants