IAP GITLAB

vertical_EAS.cc 11.1 KB
Newer Older
Ralf Ulrich's avatar
Ralf Ulrich committed
1 2 3 4 5 6 7 8 9 10 11 12
/*
 * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
 *
 * See file AUTHORS for a list of contributors.
 *
 * This software is distributed under the terms of the GNU General Public
 * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
 * the license.
 */

#include <corsika/cascade/Cascade.h>
#include <corsika/environment/Environment.h>
Ralf Ulrich's avatar
Ralf Ulrich committed
13
#include <corsika/environment/FlatExponential.h>
14 15
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMagneticFieldModel.h>
Ralf Ulrich's avatar
Ralf Ulrich committed
16
#include <corsika/environment/NuclearComposition.h>
17
#include <corsika/environment/ShowerAxis.h>
18
#include <corsika/environment/SlidingPlanarExponential.h>
19
#include <corsika/environment/UniformMagneticField.h>
Ralf Ulrich's avatar
Ralf Ulrich committed
20
#include <corsika/geometry/Plane.h>
Ralf Ulrich's avatar
Ralf Ulrich committed
21
#include <corsika/geometry/Sphere.h>
Maximilian Reininghaus's avatar
Maximilian Reininghaus committed
22 23 24 25
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/StackProcess.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/interaction_counter/InteractionCounter.h>
26
#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
Maximilian Reininghaus's avatar
Maximilian Reininghaus committed
27 28
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <corsika/process/particle_cut/ParticleCut.h>
Andre Schmidt's avatar
Andre Schmidt committed
29
#include <corsika/process/track_writer/TrackWriter.h>
Maximilian Reininghaus's avatar
Maximilian Reininghaus committed
30
#include <corsika/process/pythia/Decay.h>
31
#include <corsika/process/sibyll/Decay.h>
Ralf Ulrich's avatar
Ralf Ulrich committed
32 33
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
Maximilian Reininghaus's avatar
Maximilian Reininghaus committed
34 35
#include <corsika/process/switch_process/SwitchProcess.h>
#include <corsika/process/tracking_line/TrackingLine.h>
Ralf Ulrich's avatar
Ralf Ulrich committed
36
#include <corsika/process/urqmd/UrQMD.h>
Ralf Ulrich's avatar
Ralf Ulrich committed
37
#include <corsika/random/RNGManager.h>
Maximilian Reininghaus's avatar
Maximilian Reininghaus committed
38 39 40
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
Ralf Ulrich's avatar
Ralf Ulrich committed
41 42
#include <corsika/utl/CorsikaFenv.h>

Ralf Ulrich's avatar
Ralf Ulrich committed
43
#include <iomanip>
Ralf Ulrich's avatar
Ralf Ulrich committed
44 45
#include <iostream>
#include <limits>
46
#include <string>
Ralf Ulrich's avatar
Ralf Ulrich committed
47 48 49 50 51 52 53 54 55 56 57 58 59 60
#include <typeinfo>

using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;

using namespace std;
using namespace corsika::units::si;

Ralf Ulrich's avatar
Ralf Ulrich committed
61 62
void registerRandomStreams() {
  random::RNGManager::GetInstance().RegisterRandomStream("cascade");
63
  random::RNGManager::GetInstance().RegisterRandomStream("qgran");
Ralf Ulrich's avatar
Ralf Ulrich committed
64
  random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
65
  random::RNGManager::GetInstance().RegisterRandomStream("pythia");
Ralf Ulrich's avatar
Ralf Ulrich committed
66 67 68 69 70
  random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");

  random::RNGManager::GetInstance().SeedAll();
}

71 72
int main(int argc, char** argv) {
  if (argc != 4) {
Maximilian Reininghaus's avatar
Maximilian Reininghaus committed
73
    std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV>" << std::endl;
74 75
    return 1;
  }
Ralf Ulrich's avatar
Ralf Ulrich committed
76 77
  feenableexcept(FE_INVALID);
  // initialize random number sequence(s)
Ralf Ulrich's avatar
Ralf Ulrich committed
78
  registerRandomStreams();
Ralf Ulrich's avatar
Ralf Ulrich committed
79 80

  // setup environment, geometry
81 82
  using IEnvModel = IMagneticFieldModel<IMediumModel>;
  using EnvType = Environment<IEnvModel>;
83
  EnvType env;
Ralf Ulrich's avatar
Ralf Ulrich committed
84
  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
85
  Point const center{rootCS, 0_m, 0_m, 0_m};
86 87 88 89
  //  environment::LayeredSphericalAtmosphereBuilder builder{center};
  NuclearComposition const composition{
      {particles::Code::Nitrogen, particles::Code::Oxygen},
      {0.7847f, 1.f - 0.7847f}}; // values taken from AIRES manual, Ar removed for no
90

91
  Vector const B{rootCS, 0_T, 50_uT, 0_T};
92

93 94 95 96 97 98 99
  auto const seaLevel = 6371.0_km;

  auto node1 = std::make_unique<VolumeTreeNode<IEnvModel>>(
      std::make_unique<geometry::Sphere>(center, seaLevel + 112.8_km));
  node1->SetModelProperties<UniformMagneticField<HomogeneousMedium<IEnvModel>>>(
      B, 1_g / (1_cm * 1_cm * 1e9_cm), composition);

100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
  auto node2 = std::make_unique<VolumeTreeNode<IEnvModel>>(
      std::make_unique<geometry::Sphere>(center, seaLevel + 100.0_km));
  node2->SetModelProperties<UniformMagneticField<SlidingPlanarExponential<IEnvModel>>>(
      B, center, 540.1778_g / (1_cm * 1_cm) / 772170.16_cm, -772170.16_cm, composition,
      seaLevel);

  auto node3 = std::make_unique<VolumeTreeNode<IEnvModel>>(
      std::make_unique<geometry::Sphere>(center, seaLevel + 40.0_km));
  node3->SetModelProperties<UniformMagneticField<SlidingPlanarExponential<IEnvModel>>>(
      B, center, 1305.5948_g / (1_cm * 1_cm) / 636143.04_cm, -636143.04_cm, composition,
      seaLevel);

  auto node4 = std::make_unique<VolumeTreeNode<IEnvModel>>(
      std::make_unique<geometry::Sphere>(center, seaLevel + 10.0_km));
  node4->SetModelProperties<UniformMagneticField<SlidingPlanarExponential<IEnvModel>>>(
      B, center, 1144.9069_g / (1_cm * 1_cm) / 878153.55_cm, -878153.55_cm, composition,
      seaLevel);

  auto node5 = std::make_unique<VolumeTreeNode<IEnvModel>>(
      std::make_unique<geometry::Sphere>(center, seaLevel + 4.0_km));
  node5->SetModelProperties<UniformMagneticField<SlidingPlanarExponential<IEnvModel>>>(
      B, center, 1222.6562_g / (1_cm * 1_cm) / 994186.38_cm, -994186.38_cm, composition,
      seaLevel);

  node4->AddChild(std::move(node5));
  node3->AddChild(std::move(node4));
  node2->AddChild(std::move(node3));
Andre Schmidt's avatar
Andre Schmidt committed
127
  node1->AddChild(std::move(node2));
128
  env.GetUniverse()->AddChild(std::move(node1));
Ralf Ulrich's avatar
Ralf Ulrich committed
129 130 131 132

  // setup particle stack, and add primary particle
  setup::Stack stack;
  stack.Clear();
133 134 135 136 137
  const Code beamCode = Code::Nucleus;
  unsigned short const A = std::stoi(std::string(argv[1]));
  unsigned short Z = std::stoi(std::string(argv[2]));
  auto const mass = particles::GetNucleusMass(A, Z);
  const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
138
  double theta = 0.;
139
  auto const thetaRad = theta / 180. * M_PI;
Ralf Ulrich's avatar
Ralf Ulrich committed
140 141 142 143 144

  auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
    return sqrt((Elab - m) * (Elab + m));
  };
  HEPMomentumType P0 = elab2plab(E0, mass);
145 146
  auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
    return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
Ralf Ulrich's avatar
Ralf Ulrich committed
147
  };
148

149
  auto const [px, py, pz] = momentumComponents(thetaRad, P0);
Ralf Ulrich's avatar
Ralf Ulrich committed
150 151
  auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
  cout << "input particle: " << beamCode << endl;
152
  cout << "input angles: theta=" << theta << endl;
153 154 155
  cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
       << endl;

156 157
  auto const observationHeight = 1.4_km + seaLevel;
  auto const injectionHeight = 112.75_km + seaLevel;
158 159 160
  auto const t = -observationHeight * cos(thetaRad) +
                 sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
                      si::detail::static_pow<2>(injectionHeight));
161
  Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
162 163 164
  Point const injectionPos =
      showerCore +
      Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
165 166

  std::cout << "point of injection: " << injectionPos.GetCoordinates() << std::endl;
Ralf Ulrich's avatar
Ralf Ulrich committed
167

168 169 170 171 172 173 174
  if (A != 1) {
    stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
                                 corsika::stack::MomentumVector, geometry::Point,
                                 units::si::TimeType, unsigned short, unsigned short>{
        beamCode, E0, plab, injectionPos, 0_ns, A, Z});

  } else {
175 176 177 178 179 180 181 182 183 184 185 186 187 188
    if (Z == 1) {
      stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
                                   corsika::stack::MomentumVector, geometry::Point,
                                   units::si::TimeType>{particles::Code::Proton, E0, plab,
                                                        injectionPos, 0_ns});
    } else if (Z == 0) {
      stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
                                   corsika::stack::MomentumVector, geometry::Point,
                                   units::si::TimeType>{particles::Code::Neutron, E0,
                                                        plab, injectionPos, 0_ns});
    } else {
      std::cerr << "illegal parameters" << std::endl;
      return EXIT_FAILURE;
    }
189
  }
Ralf Ulrich's avatar
Ralf Ulrich committed
190

191 192 193 194 195
  std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.02
            << std::endl;

  environment::ShowerAxis const showerAxis{injectionPos,
                                           (showerCore - injectionPos) * 1.02, env};
Ralf Ulrich's avatar
Ralf Ulrich committed
196

Ralf Ulrich's avatar
Ralf Ulrich committed
197 198 199
  // setup processes, decays and interactions

  process::sibyll::Interaction sibyll;
200
  process::interaction_counter::InteractionCounter sibyllCounted(sibyll);
201

Ralf Ulrich's avatar
Ralf Ulrich committed
202
  process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
203
  process::interaction_counter::InteractionCounter sibyllNucCounted(sibyllNuc);
204 205 206 207

  process::pythia::Decay decayPythia;

  // use sibyll decay routine for decays of particles unknown to pythia
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227
  process::sibyll::Decay decaySibyll{{
      Code::N1440Plus,
      Code::N1440MinusBar,
      Code::N1440_0,
      Code::N1440_0Bar,
      Code::N1710Plus,
      Code::N1710MinusBar,
      Code::N1710_0,
      Code::N1710_0Bar,

      Code::Pi1300Plus,
      Code::Pi1300Minus,
      Code::Pi1300_0,

      Code::KStar0_1430_0,
      Code::KStar0_1430_0Bar,
      Code::KStar0_1430_Plus,
      Code::KStar0_1430_MinusBar,
  }};

228
  decaySibyll.PrintDecayConfig();
Ralf Ulrich's avatar
Ralf Ulrich committed
229

Andre Schmidt's avatar
Test  
Andre Schmidt committed
230
  process::particle_cut::ParticleCut cut{50_GeV};
Andre Schmidt's avatar
Andre Schmidt committed
231
  process::track_writer::TrackWriter trackWriter("tracks.dat");
Ralf Ulrich's avatar
Ralf Ulrich committed
232
  process::energy_loss::EnergyLoss eLoss(showerAxis);
233
  process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Ralf Ulrich's avatar
Ralf Ulrich committed
234

235
  Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
236 237
  process::observation_plane::ObservationPlane observationLevel(
      obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
Ralf Ulrich's avatar
Ralf Ulrich committed
238 239

  // assemble all processes into an ordered process list
Ralf Ulrich's avatar
Ralf Ulrich committed
240 241

  process::UrQMD::UrQMD urqmd;
242
  process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
Ralf Ulrich's avatar
Ralf Ulrich committed
243

244
  auto sibyllSequence = sibyllNucCounted << sibyllCounted;
245 246
  process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence,
                                                       55_GeV);
247
  auto decaySequence = decayPythia << decaySibyll;
248
  auto sequence = switchProcess << decaySequence << longprof << eLoss << cut
Andre Schmidt's avatar
Andre Schmidt committed
249
                                << observationLevel << trackWriter;
Ralf Ulrich's avatar
Ralf Ulrich committed
250

Ralf Ulrich's avatar
Ralf Ulrich committed
251
  // define air shower object, run simulation
Ralf Ulrich's avatar
Ralf Ulrich committed
252
  tracking_line::TrackingLine tracking;
Ralf Ulrich's avatar
Ralf Ulrich committed
253 254
  cascade::Cascade EAS(env, tracking, sequence, stack);
  EAS.Init();
255 256

  // to fix the point of first interaction, uncomment the following two lines:
257 258
  //  EAS.SetNodes();
  //  EAS.forceInteraction();
259

Ralf Ulrich's avatar
Ralf Ulrich committed
260 261 262 263 264 265 266 267 268 269 270
  EAS.Run();

  eLoss.PrintProfile(); // print longitudinal profile

  cut.ShowResults();
  const HEPEnergyType Efinal =
      cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
  cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
       << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
  cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
       << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
271

272 273
  auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
                     urqmdCounted.GetHistogram();
Maximilian Reininghaus's avatar
Maximilian Reininghaus committed
274 275
  hists.saveLab("inthist_lab.txt");
  hists.saveCMS("inthist_cms.txt");
276

277 278
  longprof.save("longprof.txt");

Ralf Ulrich's avatar
Ralf Ulrich committed
279 280
  std::ofstream finish("finished");
  finish << "run completed without error" << std::endl;
Ralf Ulrich's avatar
Ralf Ulrich committed
281
}