Skip to content

Commit

Permalink
Merge pull request #15 from Jean1995/master
Browse files Browse the repository at this point in the history
Include pairproduction of muons as a new process
  • Loading branch information
sudojan authored Mar 19, 2019
2 parents 96875cc + d7c8955 commit f94ff97
Show file tree
Hide file tree
Showing 24 changed files with 2,337 additions and 34 deletions.
2 changes: 1 addition & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ before_script:
- make VERBOSE=1

script:
- travis_wait 60 ctest -j2 --verbose -E "(Brems|Photo|Epair)"
- travis_wait 60 ctest -j2 --verbose -E "(Brems|Photo|Epair|Mupair)"

notifications:
email: false
Expand Down
3 changes: 3 additions & 0 deletions cmake/Standalone_CMakeLists.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -248,6 +248,7 @@ IF(ADD_TESTS)
ADD_EXECUTABLE(UnitTest_Interpolant tests/Interpolant_TEST.cxx)
ADD_EXECUTABLE(UnitTest_Bremsstrahlung tests/Bremsstrahlung_TEST.cxx)
ADD_EXECUTABLE(UnitTest_Epairproduction tests/Epairproduction_TEST.cxx)
ADD_EXECUTABLE(UnitTest_Mupairproduction tests/Mupairproduction_TEST.cxx)
ADD_EXECUTABLE(UnitTest_Ionization tests/Ionization_TEST.cxx)
ADD_EXECUTABLE(UnitTest_Medium tests/Medium_TEST.cxx)
ADD_EXECUTABLE(UnitTest_Particle tests/Particle_TEST.cxx)
Expand All @@ -268,6 +269,7 @@ IF(ADD_TESTS)
TARGET_LINK_LIBRARIES(UnitTest_Ionization PROPOSAL ${gtest_LIBRARIES})
TARGET_LINK_LIBRARIES(UnitTest_Bremsstrahlung PROPOSAL ${gtest_LIBRARIES})
TARGET_LINK_LIBRARIES(UnitTest_Epairproduction PROPOSAL ${gtest_LIBRARIES})
TARGET_LINK_LIBRARIES(UnitTest_Mupairproduction PROPOSAL ${gtest_LIBRARIES})
TARGET_LINK_LIBRARIES(UnitTest_Photonuclear PROPOSAL ${gtest_LIBRARIES})
TARGET_LINK_LIBRARIES(UnitTest_Medium PROPOSAL ${gtest_LIBRARIES})
TARGET_LINK_LIBRARIES(UnitTest_Particle PROPOSAL ${gtest_LIBRARIES})
Expand All @@ -293,6 +295,7 @@ IF(ADD_TESTS)
ADD_TEST(UnitTest_EnergyCutSettings bin/UnitTest_EnergyCutSettings)
ADD_TEST(UnitTest_Interpolant bin/UnitTest_Interpolant)
ADD_TEST(UnitTest_Epairproduction bin/UnitTest_Epairproduction)
ADD_TEST(UnitTest_Mupairproduction bin/UnitTest_Epairproduction)
ADD_TEST(UnitTest_Ionization bin/UnitTest_Ionization)
ADD_TEST(UnitTest_Bremsstrahlung bin/UnitTest_Bremsstrahlung)
ADD_TEST(UnitTest_Photonuclear bin/UnitTest_Photonuclear)
Expand Down
66 changes: 66 additions & 0 deletions private/PROPOSAL/Propagator.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -1297,6 +1297,22 @@ Sector::Definition Propagator::CreateSectorDefinition(const std::string& json_ob
log_debug("No given epair_multiplier option given. Use default (%f)", sec_def_global.utility_def.epair_def.multiplier);
}

if (json_global.find("mupair_multiplier") != json_global.end())
{
if (json_global["mupair_multiplier"].is_number())
{
sec_def_global.utility_def.mupair_def.multiplier = json_global["mupair_multiplier"].get<double>();
}
else
{
log_fatal("The given mupair_multiplier option is not a double.");
}
}
else
{
log_debug("No given mupair_multiplier option given. Use default (%f)", sec_def_global.utility_def.mupair_def.multiplier);
}

if (json_global.find("lpm") != json_global.end())
{
if (json_global["lpm"].is_boolean())
Expand All @@ -1314,6 +1330,38 @@ Sector::Definition Propagator::CreateSectorDefinition(const std::string& json_ob
log_debug("No given lpm option given. Use default (true)");
}

if (json_global.find("mupair_enable") != json_global.end())
{
if (json_global["mupair_enable"].is_boolean())
{
sec_def_global.utility_def.mupair_def.mupair_enable = json_global["mupair_enable"].get<bool>();
}
else
{
log_fatal("The given mupair_enable option is not a bool.");
}
}
else
{
log_debug("No given mupair_enable option given. Use default (false)");
}

if (json_global.find("mupair_particle_output") != json_global.end())
{
if (json_global["mupair_particle_output"].is_boolean())
{
sec_def_global.utility_def.mupair_def.particle_output = json_global["mupair_particle_output"].get<bool>();
}
else
{
log_fatal("The given mupair_particle_output option is not a bool.");
}
}
else
{
log_debug("No given mupair_particle_output option given. Use default (true)");
}

if (json_global.find("exact_time") != json_global.end())
{
if (json_global["exact_time"].is_boolean())
Expand Down Expand Up @@ -1452,6 +1500,24 @@ Sector::Definition Propagator::CreateSectorDefinition(const std::string& json_ob
PhotonuclearFactory::Get().GetStringFromEnum(sec_def_global.utility_def.photo_def.parametrization).c_str());
}

if (json_global.find("mupair") != json_global.end())
{
if (json_global["mupair"].is_string())
{
std::string mupair = json_global["mupair"].get<std::string>();
sec_def_global.utility_def.mupair_def.parametrization = MupairProductionFactory::Get().GetEnumFromString(mupair);
}
else
{
log_fatal("The given mupair option is not a string.");
}
}
else
{
log_debug("The mupair option is not set. Use default %s",
MupairProductionFactory::Get().GetStringFromEnum(sec_def_global.utility_def.mupair_def.parametrization).c_str());
}

if (json_global.find("photo_shadow") != json_global.end())
{
if (json_global["photo_shadow"].is_string())
Expand Down
3 changes: 3 additions & 0 deletions private/PROPOSAL/crossection/CrossSection.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,9 @@ std::ostream& PROPOSAL::operator<<(std::ostream& os, CrossSection const& cross)
case DynamicData::NuclInt:
name = "PhotoNuclear";
break;
case DynamicData::MuPair:
name = "MuPairProduction";
break;
default:
break;
}
Expand Down
66 changes: 66 additions & 0 deletions private/PROPOSAL/crossection/MupairIntegral.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@

#include <functional>

#include "PROPOSAL/Constants.h"
#include "PROPOSAL/crossection/MupairIntegral.h"
#include "PROPOSAL/crossection/parametrization/MupairProduction.h"
#include "PROPOSAL/medium/Medium.h"

using namespace PROPOSAL;

MupairIntegral::MupairIntegral(const MupairProduction& param)
: CrossSectionIntegral(DynamicData::MuPair, param)
{
}

MupairIntegral::MupairIntegral(const MupairIntegral& mupair)
: CrossSectionIntegral(mupair)
{
}

MupairIntegral::~MupairIntegral() {}

// ----------------------------------------------------------------- //
// Public methods
// ----------------------------------------------------------------- //

double MupairIntegral::CalculatedEdx(double energy)
{
if (parametrization_->GetMultiplier() <= 0)
{
return 0;
}

double sum = 0;

for (int i = 0; i < parametrization_->GetMedium().GetNumComponents(); i++)
{
parametrization_->SetCurrentComponent(i);
Parametrization::IntegralLimits limits = parametrization_->GetIntegralLimits(energy);

sum += dedx_integral_.Integrate(
limits.vMin,
limits.vUp,
std::bind(&Parametrization::FunctionToDEdxIntegral, parametrization_, energy, std::placeholders::_1),
4);
}

return energy * sum;
}

std::vector<Particle*> MupairIntegral::CalculateProducedParticles(double energy, double energy_loss, double rnd1, double rnd2){

//Create MuPair particles
std::vector<Particle*> mupair;
mupair.push_back(new Particle(MuMinusDef::Get()));
mupair.push_back(new Particle(MuPlusDef::Get()));

//Sample and assign energies
double rho = parametrization_->Calculaterho(energy, energy_loss/energy, rnd1, rnd2);

mupair[0]->SetEnergy(0.5*energy_loss*(1 + rho));
mupair[1]->SetEnergy(0.5*energy_loss*(1 - rho));
return mupair;

}

111 changes: 111 additions & 0 deletions private/PROPOSAL/crossection/MupairInterpolant.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@

#include <functional>

#include "PROPOSAL/crossection/MupairIntegral.h"
#include "PROPOSAL/crossection/MupairInterpolant.h"
#include "PROPOSAL/crossection/parametrization/MupairProduction.h"

#include "PROPOSAL/math/Interpolant.h"
#include "PROPOSAL/math/InterpolantBuilder.h"

#include "PROPOSAL/Constants.h"
#include "PROPOSAL/Output.h"
#include "PROPOSAL/methods.h"

using namespace PROPOSAL;

MupairInterpolant::MupairInterpolant(const MupairProduction& param, InterpolationDef def)
: CrossSectionInterpolant(DynamicData::MuPair, param)
{
// Use parent CrossSecition dNdx interpolation
InitdNdxInerpolation(def);

// --------------------------------------------------------------------- //
// Builder for DEdx
// --------------------------------------------------------------------- //

Interpolant1DBuilder builder1d;
Helper::InterpolantBuilderContainer builder_container;

// Needed for CalculatedEdx integration
MupairIntegral mupair(param);

builder1d.SetMax(NUM1)
.SetXMin(param.GetParticleDef().low)
.SetXMax(BIGENERGY)
.SetRomberg(def.order_of_interpolation)
.SetRational(true)
.SetRelative(false)
.SetIsLog(true)
.SetRombergY(def.order_of_interpolation)
.SetRationalY(false)
.SetRelativeY(false)
.SetLogSubst(true)
.SetFunction1D(std::bind(&CrossSection::CalculatedEdx, &mupair, std::placeholders::_1));

builder_container.push_back(std::make_pair(&builder1d, &dedx_interpolant_));

// --------------------------------------------------------------------- //
// Builder for DE2dx
// --------------------------------------------------------------------- //

Interpolant1DBuilder builder_de2dx;
Helper::InterpolantBuilderContainer builder_container_de2dx;

builder_de2dx.SetMax(NUM2)
.SetXMin(param.GetParticleDef().low)
.SetXMax(BIGENERGY)
.SetRomberg(def.order_of_interpolation)
.SetRational(false)
.SetRelative(false)
.SetIsLog(true)
.SetRombergY(def.order_of_interpolation)
.SetRationalY(false)
.SetRelativeY(false)
.SetLogSubst(false)
.SetFunction1D(std::bind(&CrossSection::CalculatedE2dx, &mupair, std::placeholders::_1));

builder_container_de2dx.push_back(std::make_pair(&builder_de2dx, &de2dx_interpolant_));

Helper::InitializeInterpolation("dEdx", builder_container, std::vector<Parametrization*>(1, parametrization_), def);
Helper::InitializeInterpolation(
"dE2dx", builder_container_de2dx, std::vector<Parametrization*>(1, parametrization_), def);
}

MupairInterpolant::MupairInterpolant(const MupairInterpolant& mupair)
: CrossSectionInterpolant(mupair)
{
}

MupairInterpolant::~MupairInterpolant() {}

// ----------------------------------------------------------------- //
// Public methods
// ----------------------------------------------------------------- //

// ------------------------------------------------------------------------- //
double MupairInterpolant::CalculatedEdx(double energy)
{
if (parametrization_->GetMultiplier() <= 0)
{
return 0;
}

return std::max(dedx_interpolant_->Interpolate(energy), 0.0);
}

std::vector<Particle*> MupairInterpolant::CalculateProducedParticles(double energy, double energy_loss, double rnd1, double rnd2){

//Create MuPair particles
std::vector<Particle*> mupair;
mupair.push_back(new Particle(MuMinusDef::Get()));
mupair.push_back(new Particle(MuPlusDef::Get()));

//Sample and assign energies
double rho = parametrization_->Calculaterho(energy, energy_loss/energy, rnd1, rnd2);

mupair[0]->SetEnergy(0.5*energy_loss*(1 + rho));
mupair[1]->SetEnergy(0.5*energy_loss*(1 - rho));
return mupair;

}
Loading

0 comments on commit f94ff97

Please sign in to comment.