From da515ab451c62cbce4b8947ab1ca15a0d7791e62 Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Tue, 20 Apr 2021 18:18:49 +0200
Subject: [PATCH] Synchrotron radiation using C8 set up
---
examples/radio_shower2.cpp | 97 ++++++++++++++++++----------
tests/modules/testRadio.cpp | 122 ++++++++++++++++++------------------
2 files changed, 125 insertions(+), 94 deletions(-)
diff --git a/examples/radio_shower2.cpp b/examples/radio_shower2.cpp
index 3f8d5353b..ad0440ae5 100644
--- a/examples/radio_shower2.cpp
+++ b/examples/radio_shower2.cpp
@@ -40,6 +40,7 @@
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/TimeCut.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/HadronicElasticModel.hpp>
#include <corsika/modules/Pythia8.hpp>
@@ -70,73 +71,104 @@ int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
- std::cout << "radio_shower2" << std::endl;
+ std::cout << "Synchrotron radiation" << std::endl;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
RNGManager::getInstance().registerRandomStream("cascade");
- // This environment needs a hardcoded refractive index in the propagator at the moment
- // setup environment, geometry
+ // This environment may need a hardcoded refractive index in the propagator at the moment although it shouldn't
+ // set up the environment
using EnvType = setup::Environment;
EnvType env;
auto& universe = *(env.getUniverse());
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
- // the antenna location
- const auto point1{Point(rootCS, 50_m, 50_m, 0_m)};
- const auto point2{Point(rootCS, 50_m, -50_m, 0_m)};
- const auto point3{Point(rootCS, -50_m, 50_m, 0_m)};
- const auto point4{Point(rootCS, -50_m, -50_m, 0_m)};
-
- // the antennas
- TimeDomainAntenna ant1("antenna1", point1, 0_s, 1_s, 1/1e-6_s);
- TimeDomainAntenna ant2("antenna2", point2, 0_s, 1_s, 1/1e-6_s);
- TimeDomainAntenna ant3("antenna3", point3, 0_s, 1_s, 1/1e-6_s);
- TimeDomainAntenna ant4("antenna4", point4, 0_s, 1_s, 1/1e-6_s);
-
- // the detector
- AntennaCollection<TimeDomainAntenna> detector;
- detector.addAntenna(ant1);
- detector.addAntenna(ant2);
- detector.addAntenna(ant3);
- detector.addAntenna(ant4);
-
auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
using MyHomogeneousModel = UniformRefractiveIndex<MediumPropertyModel<
UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>>;
- world->setModelProperties<MyHomogeneousModel>(1.000327,
- Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 1_T),
+ world->setModelProperties<MyHomogeneousModel>(1,
+ Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 0.3809_T),
1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{Code::Hydrogen},
+ NuclearComposition(std::vector<Code>{Code::Nitrogen},
std::vector<float>{(float)1.}));
universe.addChild(std::move(world));
+// // The following environment is the same as the one above qualitatively but it doesn't compile with Cascade.
+// // I leave it here for now, as I am curious to understand why at some point.
+// using IModelInterface = IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+// using AtmModel = UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<HomogeneousMedium
+// <IModelInterface>>>>;
+// using EnvType = Environment<AtmModel>;
+// EnvType env;
+// auto& universe = *(env.getUniverse());
+// CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+//
+// auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
+//
+// world->setModelProperties<AtmModel>(1,
+// Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 0.3809_T),
+// 1_kg / (1_m * 1_m * 1_m),
+// NuclearComposition(std::vector<Code>{Code::Nitrogen},
+// std::vector<float>{(float)1.}));
+//
+// universe.addChild(std::move(world));
+
+ // the antenna locations
+ const auto point1{Point(rootCS, 100_m, 100_m, 0_m)};
+ const auto point2{Point(rootCS, 100_m, -100_m, 0_m)};
+ const auto point3{Point(rootCS, -100_m, -100_m, 0_m)};
+ const auto point4{Point(rootCS, -100_m, 100_m, 0_m)};
+
+ // the antenna time variables
+ const TimeType t1{0_s};
+ const TimeType t2{1e-6_s};
+ const InverseTimeType t3{1e+9_Hz};
+
+ // the antennas
+ TimeDomainAntenna ant1("antenna 1", point1, t1, t2, t3);
+ TimeDomainAntenna ant2("antenna 2", point2, t1, t2, t3);
+ TimeDomainAntenna ant3("antenna 3", point3, t1, t2, t3);
+ TimeDomainAntenna ant4("antenna 4", point4, t1, t2, t3);
+
+ // the detector
+ AntennaCollection<TimeDomainAntenna> detector;
+ detector.addAntenna(ant1);
+ detector.addAntenna(ant2);
+ detector.addAntenna(ant3);
+ detector.addAntenna(ant4);
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
const Code beamCode = Code::Electron;
const HEPMassType mass = Electron::mass;
- const HEPEnergyType E0 = 1000_GeV;
+ const HEPEnergyType E0 = 11.4_MeV;
+ double theta = 0.;
+ double phi = 0.;
- Point injectionPos(rootCS, 0_m, 0_m, 0_m);
+ Point injectionPos(rootCS, 0_m, 100_m, 0_m);
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
};
HEPMomentumType P0 = elab2plab(E0, mass);
-
- auto plab = MomentumVector(rootCS, {0, P0, 0});
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * cos(theta), ptot * sin(theta),
+ ptot * sin(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ auto plab = MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
}
- // setup processes, decays and interactions
+ // setup relevant processes
setup::Tracking tracking;
// StackInspector<setup::Stack> stackInspect(1, true, E0);
@@ -145,18 +177,17 @@ int main() {
decltype(StraightPropagator(env))>, decltype(StraightPropagator(env))>
coreas(detector, env);
+ TimeCut cut(1e-9_s);
TrackWriter trackWriter("tracks.dat");
-// ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
// assemble all processes into an ordered process list
- auto sequence = make_sequence(coreas, trackWriter);
+ auto sequence = make_sequence(coreas, cut, trackWriter);
// define air shower object, run simulation
Cascade EAS(env, tracking, sequence, stack);
EAS.run();
- //TODO: this will run indefinetly due to no energy losses
// get radio output
coreas.writeOutput();
}
diff --git a/tests/modules/testRadio.cpp b/tests/modules/testRadio.cpp
index 7d5ff6e31..7b1ba2be7 100644
--- a/tests/modules/testRadio.cpp
+++ b/tests/modules/testRadio.cpp
@@ -298,12 +298,12 @@ TEST_CASE("Radio", "[processes]") {
env.getUniverse()->addChild(std::move(Medium));
- // now create antennas and detectors/////////////////////////////////////////////
+ // now create antennas and detectors/////////////////////////////////////////////
// the antennas location
- const auto point1{Point(rootCS, 200_m, 200_m, 0_m)};
- const auto point2{Point(rootCS, 200_m, -200_m, 0_m)};
- const auto point3{Point(rootCS, -200_m, -200_m, 0_m)};
- const auto point4{Point(rootCS, -200_m, 200_m, 0_m)};
+ const auto point1{Point(rootCS, 100_m, 100_m, 0_m)};
+ const auto point2{Point(rootCS, 100_m, -100_m, 0_m)};
+ const auto point3{Point(rootCS, -100_m, -100_m, 0_m)};
+ const auto point4{Point(rootCS, -100_m, 100_m, 0_m)};
// create times for the antenna
@@ -734,45 +734,45 @@ TEST_CASE("Radio", "[processes]") {
// store all the points in a std::array
std::array<Point, 400> points_
{p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,
- p10,p11,p12,p13,p14,p15,p16,p17,p18,p19,
- p20,p21,p22,p23,p24,p25,p26,p27,p28,p29,
- p30,p31,p32,p33,p34,p35,p36,p37,p38,p39,
- p40,p41,p42,p43,p44,p45,p46,p47,p48,p49,
- p50,p51,p52,p53,p54,p55,p56,p57,p58,p59,
- p60,p61,p62,p63,p64,p65,p66,p67,p68,p69,
- p70,p71,p72,p73,p74,p75,p76,p77,p78,p79,
- p80,p81,p82,p83,p84,p85,p86,p87,p88,p89,
- p90,p91,p92,p93,p94,p95,p96,p97,p98,p99,
- p100,p101,p102,p103,p104,p105,p106,p107,p108,p109,
- p110,p111,p112,p113,p114,p115,p116,p117,p118,p119,
- p120,p121,p122,p123,p124,p125,p126,p127,p128,p129,
- p130,p131,p132,p133,p134,p135,p136,p137,p138,p139,
- p140,p141,p142,p143,p144,p145,p146,p147,p148,p149,
- p150,p151,p152,p153,p154,p155,p156,p157,p158,p159,
- p160,p161,p162,p163,p164,p165,p166,p167,p168,p169,
- p170,p171,p172,p173,p174,p175,p176,p177,p178,p179,
- p180,p181,p182,p183,p184,p185,p186,p187,p188,p189,
- p190,p191,p192,p193,p194,p195,p196,p197,p198,p199,
- p200,p201,p202,p203,p204,p205,p206,p207,p208,p209,
- p210,p211,p212,p213,p214,p215,p216,p217,p218,p219,
- p220,p221,p222,p223,p224,p225,p226,p227,p228,p229,
- p230,p231,p232,p233,p234,p235,p236,p237,p238,p239,
- p240,p241,p242,p243,p244,p245,p246,p247,p248,p249,
- p250,p251,p252,p253,p254,p255,p256,p257,p258,p259,
- p260,p261,p262,p263,p264,p265,p266,p267,p268,p269,
- p270,p271,p272,p273,p274,p275,p276,p277,p278,p279,
- p280,p281,p282,p283,p284,p285,p286,p287,p288,p289,
- p290,p291,p292,p293,p294,p295,p296,p297,p298,p299,
- p300,p301,p302,p303,p304,p305,p306,p307,p308,p309,
- p310,p311,p312,p313,p314,p315,p316,p317,p318,p319,
- p320,p321,p322,p323,p324,p325,p326,p327,p328,p329,
- p330,p331,p332,p333,p334,p335,p336,p337,p338,p339,
- p340,p341,p342,p343,p344,p345,p346,p347,p348,p349,
- p350,p351,p352,p353,p354,p355,p356,p357,p358,p359,
- p360,p361,p362,p363,p364,p365,p366,p367,p368,p369,
- p370,p371,p372,p373,p374,p375,p376,p377,p378,p379,
- p380,p381,p382,p383,p384,p385,p386,p387,p388,p389,
- p390,p391,p392,p393,p394,p395,p396,p397,p398,p399};
+ p10,p11,p12,p13,p14,p15,p16,p17,p18,p19,
+ p20,p21,p22,p23,p24,p25,p26,p27,p28,p29,
+ p30,p31,p32,p33,p34,p35,p36,p37,p38,p39,
+ p40,p41,p42,p43,p44,p45,p46,p47,p48,p49,
+ p50,p51,p52,p53,p54,p55,p56,p57,p58,p59,
+ p60,p61,p62,p63,p64,p65,p66,p67,p68,p69,
+ p70,p71,p72,p73,p74,p75,p76,p77,p78,p79,
+ p80,p81,p82,p83,p84,p85,p86,p87,p88,p89,
+ p90,p91,p92,p93,p94,p95,p96,p97,p98,p99,
+ p100,p101,p102,p103,p104,p105,p106,p107,p108,p109,
+ p110,p111,p112,p113,p114,p115,p116,p117,p118,p119,
+ p120,p121,p122,p123,p124,p125,p126,p127,p128,p129,
+ p130,p131,p132,p133,p134,p135,p136,p137,p138,p139,
+ p140,p141,p142,p143,p144,p145,p146,p147,p148,p149,
+ p150,p151,p152,p153,p154,p155,p156,p157,p158,p159,
+ p160,p161,p162,p163,p164,p165,p166,p167,p168,p169,
+ p170,p171,p172,p173,p174,p175,p176,p177,p178,p179,
+ p180,p181,p182,p183,p184,p185,p186,p187,p188,p189,
+ p190,p191,p192,p193,p194,p195,p196,p197,p198,p199,
+ p200,p201,p202,p203,p204,p205,p206,p207,p208,p209,
+ p210,p211,p212,p213,p214,p215,p216,p217,p218,p219,
+ p220,p221,p222,p223,p224,p225,p226,p227,p228,p229,
+ p230,p231,p232,p233,p234,p235,p236,p237,p238,p239,
+ p240,p241,p242,p243,p244,p245,p246,p247,p248,p249,
+ p250,p251,p252,p253,p254,p255,p256,p257,p258,p259,
+ p260,p261,p262,p263,p264,p265,p266,p267,p268,p269,
+ p270,p271,p272,p273,p274,p275,p276,p277,p278,p279,
+ p280,p281,p282,p283,p284,p285,p286,p287,p288,p289,
+ p290,p291,p292,p293,p294,p295,p296,p297,p298,p299,
+ p300,p301,p302,p303,p304,p305,p306,p307,p308,p309,
+ p310,p311,p312,p313,p314,p315,p316,p317,p318,p319,
+ p320,p321,p322,p323,p324,p325,p326,p327,p328,p329,
+ p330,p331,p332,p333,p334,p335,p336,p337,p338,p339,
+ p340,p341,p342,p343,p344,p345,p346,p347,p348,p349,
+ p350,p351,p352,p353,p354,p355,p356,p357,p358,p359,
+ p360,p361,p362,p363,p364,p365,p366,p367,p368,p369,
+ p370,p371,p372,p373,p374,p375,p376,p377,p378,p379,
+ p380,p381,p382,p383,p384,p385,p386,p387,p388,p389,
+ p390,p391,p392,p393,p394,p395,p396,p397,p398,p399};
std::vector<TimeType> times_;
@@ -803,11 +803,11 @@ TEST_CASE("Radio", "[processes]") {
StraightTrajectory track {l,t};
auto particle1{stack.addParticle(std::make_tuple(particle, E0, plab, points_[i-1], t))}; //TODO: plab is inconsistent
coreas.doContinuous(particle1,track,true);
- }
+ }
- // ToDo: use just one electron, this way I have 400! (although it shouldn't affect the physics)
+ // ToDo: use just one electron, this way I have 400! (although it shouldn't affect the physics)
- // get the last track
+ // get the last track
TimeType t {(points_[0] - points_[399]).getNorm() / (0.999 * constants::c)};
VelocityVector v { (points_[0] - points_[399]) / t };
auto beta {v / constants::c};
@@ -819,7 +819,7 @@ TEST_CASE("Radio", "[processes]") {
coreas.doContinuous(particle1,track,true);
// get the output
- coreas.writeOutput();
+ coreas.writeOutput();
}
@@ -943,12 +943,12 @@ TEST_CASE("Radio", "[processes]") {
}
- // check that I can create working Straight Propagators in different environments
+ // check that I can create working Straight Propagators in different environments
SECTION("Straight Propagator w/ Uniform Refractive Index") {
// create an environment with uniform refractive index of 1
using UniRIndex =
- UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
+ UniformRefractiveIndex<HomogeneousMedium<IRefractiveIndexModel<IMediumModel>>>;
using EnvType = Environment<IRefractiveIndexModel<IMediumModel>>;
EnvType env;
@@ -1018,12 +1018,12 @@ TEST_CASE("Radio", "[processes]") {
CHECK(paths_.size() == 1);
}
- SECTION("Straight Propagator w/ Exponential Refractive Index") {
+ SECTION("Straight Propagator w/ Exponential Refractive Index") {
// logging::set_level(logging::level::info);
// corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
- // create an environment with exponential refractive index (n_0 = 1 & lambda = 0)
+ // create an environment with exponential refractive index (n_0 = 1 & lambda = 0)
using ExpoRIndex = ExponentialRefractiveIndex<HomogeneousMedium
<IRefractiveIndexModel<IMediumModel>>>;
@@ -1039,10 +1039,10 @@ TEST_CASE("Radio", "[processes]") {
auto const props1 =
Medium1
->setModelProperties<ExpoRIndex>( 1, 0 / 1_m,
- 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(
- std::vector<Code>{Code::Nitrogen},
- std::vector<float>{1.f}));
+ 1_kg / (1_m * 1_m * 1_m),
+ NuclearComposition(
+ std::vector<Code>{Code::Nitrogen},
+ std::vector<float>{1.f}));
env1.getUniverse()->addChild(std::move(Medium1));
@@ -1137,9 +1137,9 @@ TEST_CASE("Radio", "[processes]") {
CHECK( path.receive_.getComponents() == vvv2.getComponents() );
CHECK( path.R_distance_ == 10_m );
}
+//
+// CHECK( paths2_.size() == 1 );
+//
+// }
- CHECK( paths2_.size() == 1 );
-
- }
-
-} // END: TEST_CASE("Radio", "[processes]")
+ } // END: TEST_CASE("Radio", "[processes]")
--
GitLab