Muon Stopping Power and Range

Each example consists of a main code to be linked with the Muon Range Library. The examples calculate the range or energy loss in materials, the material thickness that would cause a given energy decrement, or simply output a formatted stopping power and range table.

Example Codes

After package installation, the example codes can be found in /usr/share/doc/murange/examples/.

Table of Contents

Example 1: Calculate the range of a 200 MeV/c muon and 1 GeV kinetic energy in Iron
Example 2: Print a stopping power and range table for Stainless Steel 316L
Example 3: Calculate the energy loss of a 30 GeV muon in 5 g/cm² Deuterium gas at 5 atm
Example 4: Calculate the energy loss of a 10 MeV muon in 1 g/cm² Cobalt-Chromium allow (ASTM-F75)
Example 5: Calculate the thickness of concrete needed to reduce the muon energy from 300 MeV to 200 MeV
Example 6: Calculate energy range of a 30 GeV muon in Carbon Tetrafluoride at 1 atm and 100 Torr
Example 7: Calculate energy loss iteratively to estimate the computing performance

← Return to Muon Range Library

Example 1

The following example code creates an element, a range table for the element, a muon object, and outputs the range at muon momentum 200 MeV/c and 1 GeV kinetic energy,

using namespace std;

#include <iostream>
#include <murange.hpp>

int main() {

  // Iron (Fe)
  Material *Fe = new Element("Fe");

  // range table for element Fe
  Range *rFe = new Range(Fe);

  // muon with momentum 200 MeV/c
  Muon *mu = new Muon();
  mu->SetMomentum(200);

  // method 1
  cout << rFe->GetRange(mu) << " MeV/cm2" << endl;

  // method 2
  cout << rFe->GetRange(1000) << " MeV/cm2" << endl;

  delete mu;
  delete rFe;
  delete Fe;

  return EXIT_SUCCESS;
}

The outcome of this example is,

$ g++ -o GetRange-example GetRange-example.cpp -lmurange
$ ./GetRange-example
56.7843 MeV/cm2
640.323 MeV/cm2

← Return to Muon Range Library

Example 2

The following example code creates a material (Stainless Steel 316L) and prints to standard output a stopping power and range table,

using namespace std;

#include <iostream>
#include <murange.hpp>

int main() {

  // Stainless Steel 316L
  Material *mat = new Mixture(502);
  Range *rng = new Range(mat);

  // print stopping power and range table
  rng->PrintRangeTable();

  delete rng;
  delete mat;

  return EXIT_SUCCESS;
}

The outcome of this example is,

$ g++ -o PrintRangeTable-example PrintRangeTable-example.cpp -lmurange
$ ./PrintRangeTable-example
̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶
 Stainless Steel (316L) (502)

   <Z/A>    ρ[g/cm³]     I[ev]      a       k=ms         x₀         x₁         C       δ₀
  0.46537   7.990000     323.4   0.16410   3.00000     0.2000     3.0000     4.5224   0.00
̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶
         T        p     Ioization      Brems     Pair prod   Photonucl    Total     CSDA Range
       [MeV]   [MeV/c]  [MeV cm²/g] [MeV cm²/g] [MeV cm²/g] [MeV cm²/g] [MeV cm²/g]   [g/cm²]
̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶
        10.0  4.704E+01    5.415       0.000       0.000       0.000       5.415     1.054E+00
        14.0  5.616E+01    4.261       0.000       0.000       0.000       4.261     1.895E+00
        20.0  6.802E+01    3.355       0.000       0.000       0.000       3.355     3.500E+00
        30.0  8.509E+01    2.625       0.000       0.000       0.000       2.625     6.915E+00
        40.0  1.003E+02    2.254       0.000       0.000       0.000       2.254     1.105E+01
        80.0  1.527E+02    1.714       0.000       0.000       0.000       1.714     3.203E+01

       100.0  1.764E+02    1.617       0.000       0.000       0.000       1.617     4.407E+01
       140.0  2.218E+02    1.519       0.000       0.000       0.000       1.519     6.971E+01
       200.0  2.868E+02    1.467       0.000       0.000       0.000       1.467     1.100E+02
       300.0  3.917E+02    1.457       0.000       0.000       0.000       1.457     1.786E+02
       400.0  4.945E+02    1.473       0.000       0.000       0.000       1.473     2.469E+02
       800.0  8.995E+02    1.554       0.000       0.000       0.000       1.554     5.112E+02

      1000.0  1.101E+03    1.587       0.001       0.000       0.000       1.588     6.384E+02
      1400.0  1.502E+03    1.641       0.001       0.000       0.001       1.643     8.858E+02
      2000.0  2.103E+03    1.700       0.002       0.001       0.001       1.703     1.244E+03
      3000.0  3.104E+03    1.765       0.003       0.002       0.001       1.771     1.819E+03
      4000.0  4.104E+03    1.811       0.004       0.004       0.002       1.820     2.375E+03
      8000.0  8.105E+03    1.914       0.010       0.011       0.003       1.938     4.498E+03

     10000.0  1.011E+04    1.944       0.014       0.015       0.004       1.977     5.520E+03
     14000.0  1.411E+04    1.989       0.021       0.024       0.006       2.040     7.510E+03
     20000.0  2.011E+04    2.033       0.033       0.039       0.008       2.113     1.040E+04
     30000.0  3.011E+04    2.080       0.054       0.067       0.012       2.214     1.502E+04
     40000.0  4.011E+04    2.112       0.077       0.098       0.016       2.302     1.945E+04
     80000.0  8.011E+04    2.183       0.172       0.235       0.031       2.621     3.573E+04

    100000.0  1.001E+05    2.205       0.223       0.309       0.038       2.775     4.314E+04
    140000.0  1.401E+05    2.237       0.328       0.462       0.054       3.081     5.683E+04
    200000.0  2.001E+05    2.271       0.492       0.704       0.076       3.543     7.498E+04
    300000.0  3.001E+05    2.310       0.772       1.110       0.114       4.306     1.006E+05
    400000.0  4.001E+05    2.337       1.062       1.530       0.152       5.081     1.220E+05
    800000.0  8.001E+05    2.404       2.263       3.259       0.307       8.232     1.837E+05

   1000000.0  1.000E+06    2.426       2.880       4.144       0.386       9.835     2.059E+05
   1400000.0  1.400E+06    2.459       4.119       5.912       0.547      13.036     2.412E+05
   2000000.0  2.000E+06    2.494       6.017       8.612       0.791      17.913     2.803E+05
   3000000.0  3.000E+06    2.535       9.188      13.102       1.210      26.034     3.265E+05
   4000000.0  4.000E+06    2.564      12.404      17.641       1.636      34.245     3.599E+05
   8000000.0  8.000E+06    2.636      25.386      35.893       3.404      67.319     4.421E+05

  10000000.0  1.000E+07    2.660      31.934      45.072       4.313      83.979     4.687E+05
  14000000.0  1.400E+07    2.696      44.995      63.382       6.180     117.253     5.088E+05
  20000000.0  2.000E+07    2.735      64.715      90.976       9.045     167.471     5.514E+05
  30000000.0  3.000E+07    2.780      97.544     136.895      14.000     251.220     5.999E+05
  40000000.0  4.000E+07    2.813     130.505     182.930      19.077     335.325     6.342E+05
  80000000.0  8.000E+07    2.893     262.586     367.168      40.303     672.950     7.167E+05

 100000000.0  1.000E+08    2.920     328.757     459.361      51.277     842.314     7.432E+05
̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶

← Return to Muon Range Library

Example 3

The following example code creates an element (Deuterium, gaseous), a range table for the element, a muon object with 30 GeV kinetic energy, and outputs the energy loss and outgoing energy after passage through 5 g/cm² D gas at 5 atm,

using namespace std;

#include <iostream>
#include <murange.hpp>

int main() {

  // Deuterium (D)
  Material *D = new Element("D");

  // set density to 5 atm
  D->SetDensity(5,"atm");

  // range table for gaseous D
  Range *rD = new Range(D);

  // muon with kinetic energy 30 GeV
  Muon *mu = new Muon(30000);

  // method 1
  D->SetThickness(5);
  cout << rD->GetOutgoingEnergy(mu) << " MeV energy loss" << endl;
  cout << mu->GetKEnergy() << " MeV" << endl;
  cout << endl;

  // method 2
  mu->SetKEnergy(30000);
  cout << rD->GetOutgoingEnergy(mu,5) << " MeV energy loss" << endl;
  cout << mu->GetKEnergy() << " MeV" << endl;
  cout << endl;

  // method 3
  cout << rD->GetOutgoingEnergy(30000,5) << " MeV energy loss" << endl;

  delete mu;
  delete rD;
  delete D;

  return EXIT_SUCCESS;
}

The outcome of this example is,

$ g++ -o GetOutgoingEnergy-example GetOutgoingEnergy-example.cpp -lmurange
$ ./GetOutgoingEnergy-example
14.9812 MeV energy loss
29985 MeV

14.9812 MeV energy loss
29985 MeV

14.9812 MeV energy loss

← Return to Muon Range Library

Example 4

The following example code creates a mixture (Co-Cr ASTM-F75), a range table for the mixture, a muon object with 10 MeV kinetic energy, and outputs the energy loss and ingoing energy before passage through 1 g/cm² liquid Deuterium,

using namespace std;

#include <iostream>
#include <murange.hpp>

int main() {

  // muon with kinetic energy 10 MeV
  Muon *mu = new Muon(10);

  // Cobalt-Chromium (ASTM-F75)
  Material *CoCr = new Mixture(532);

  // range table for Cobalt-Chromium
  Range *rCoCr = new Range(CoCr);

  // method 1
  CoCr->SetThickness(1);
  cout << rCoCr->GetIngoingEnergy(mu) << " MeV energy loss" << endl;
  cout << mu->GetKEnergy() << " MeV" << endl;
  cout << endl;

  // method 2
  mu->SetKEnergy(10);
  cout << rCoCr->GetIngoingEnergy(mu,1) << " MeV energy loss" << endl;
  cout << mu->GetKEnergy() << " MeV" << endl;
  cout << endl;

  // method 3
  cout << rCoCr->GetIngoingEnergy(10,1) << " MeV energy loss" << endl;

  delete rCoCr;
  delete CoCr;
  delete mu;

  return EXIT_SUCCESS;
}

The outcome of this example is,

$ g++ -o GetIngoingEnergy-example GetIngoingEnergy-example.cpp -lmurange
$ ./GetIngoingEnergy-example
4.59922 MeV energy loss
14.5992 MeV

4.59922 MeV energy loss
14.5992 MeV

4.59922 MeV energy loss

← Return to Muon Range Library

Example 5

The following example code creates an mixture (Concrete), a range table for the mixture, a muon object with 300 MeV kinetic energy, and outputs the material thickness needed to reduce the kinetic energy to 200 MeV.

using namespace std;

#include <iostream>
#include <murange.hpp>

int main() {

  // Concrete (Ordinary)
  Material *mat = new Mixture(518);

  // range table for concrete
  Range *r = new Range(mat);

  // muon with kinetic energy 300 MeV
  Muon *mu = new Muon(300.);

  // method 1
  double s = r->GetMaterialThickness(mu,-100.);
  cout << "material thickness " << s << " g/cm2, or " << s/mat->GetDensity() << " cm concrete" << endl;
  cout << "initial KE " << mu->GetKEnergy() << " MeV, energy loss " << r->GetOutgoingEnergy(mu,s) <<
    " MeV, final KE " << mu->GetKEnergy() << " MeV" << endl;

  // method 2
  s = r->GetMaterialThickness(300.,-100.);
  cout << "material thickness " << s << " g/cm2, or " << s/mat->GetDensity() << " cm concrete" << endl;

  delete mu;
  delete r;
  delete mat;

  return EXIT_SUCCESS;
}

The outcome of this example is,

$ g++ -o GetMaterialThickness-example GetMaterialThickness-example.cpp -lmurange
$ ./GetMaterialThickness-example
material thickness 57.4584 g/cm2, or 24.9819 cm concrete
initial KE 300 MeV, energy loss 100 MeV, final KE 200 MeV
material thickness 57.4584 g/cm2, or 24.9819 cm concrete

← Return to Muon Range Library

Example 6

The following example code calculates the range of a 30 GeV muon in Carbon Tetrafluoride (CF4) at 1 atm and 100 Torr,

using namespace std;

#include <iostream>
#include <murange.hpp>

int main() {

  // muon with kinetic energy 30 GeV
  Muon *mu = new Muon(30e3);

  // range table for CF4
  Material *CF4 = new Compound(326);
  Range *rCF4 = new Range(CF4);

  cout << "Range of" << mu->GetKEnergy()/1000 << " GeV muon in CF4 at 1 atm: " <<
    rCF4->GetRange(mu) << " MeV/cm2" << endl;

  // change the density to 100 Torr
  CF4->SetDensity(100,"Torr");

  // new range table for CF4
  delete rCF4;
  rCF4 = new Range(CF4);

  cout << "Range of" << mu->GetKEnergy()/1000 << " GeV muon in CF4 at 100 Torr: " <<
    rCF4->GetRange(mu) << " MeV/cm2" << endl;

  delete rCF4;
  delete CF4;
  delete mu;

  return EXIT_SUCCESS;
}

The outcome of this example is,

$ g++ -o cf4-range cf4-range.cpp -lmurange
$ ./cf4-range
Range of 30 GeV muon in CF4 at 1 atm: 12171.4 MeV/cm2
Range of 30 GeV muon in CF4 at 1 Torr: 11934.8 MeV/cm2

This demonstrates the density dependence of the electronic stopping power.

← Return to Muon Range Library

Example 7

The following example code creates an element (Si), a range table and sets the muon kinetic energy randomly, from 10 to 100 MeV, then calculates the energy loss in 1-10 g/cm² Si, iterated over tries times. The OpenMP time function is used to estimate the performance [Mtrials/s],

using namespace std;

#include <iostream>
#include <ctime>
#include <murange.hpp>

#include <omp.h>

int main(int argc, char **argv) {

  // number of tries
  unsigned long tries = atol(argv[1]);

  srand(time(NULL));

  // start time
  double start_time = omp_get_wtime();

  // Si range table
  Material *Si = new Element("Si");
  Range *rSi = new Range(Si);

  for ( int i = 0 ; i < tries ; i++ ) {
    // kinetic energy between 10 and 100 MeV
    double ke = (rand()%(90000000)+10000000)/1000000.;
    // thickness between 1 and 10 g/cm2
    double s = (rand()%(90000000)+10000000)/10000000.;
    rSi->GetOutgoingEnergy(ke,s);
  }

  // stop time
  double stop_time = omp_get_wtime();

  double performance = (double)tries / (time_stop-time_start) / 1000000.;

  fprintf(stderr,"\nnumber of tries: %ld\n",tries);
  fprintf(stderr,"execution time: %g s\n",time_stop-time_start);
  fprintf(stderr,"performance: %g Mtrials/s\n\n",performance);

  delete rSi;
  delete Si;

  return EXIT_SUCCESS;;
}

The outcome of this example is,

$ g++ -o performance performance.cpp -fopenmp -lmurange
$ ./performance 100000000

number of tries: 100000000
execution time: 18.8887 s
performance: 5.29418 Mtrials/s

The plot shows the performance as a function of number of trials, a useful measure when performing Monte Carlo simulations.

A similar OpenMP version of the above code, with 8 threads, gives,

$ g++ -o performance-omp performance-omp.cpp -fopenmp -lmurange
$ ./performance-omp 8 100000000
using OpenMP with 8/16 threads

number of tries: 100000000
execution time: 3.22482 s
performance: 31.0095 Mtrials/s

and performance plot,

← Return to Muon Range Library