Skip to content
Snippets Groups Projects
Commit 00922e02 authored by Steffo's avatar Steffo
Browse files

initial

parents
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
bethebloch.hpp 0 → 100644
+ 316
− 0
View file @ 00922e02
/* **************************************************************** */
/* */
/* */
/* Ziel : Berechnung des Energieverlustes nach Bethe - Bloch */
/* */
/* Datum : 01/2021 */
/* */
/* Autor : Stephan Bökelmann, auf Basis von Miriam Fritsch */
/* */
/* **************************************************************** */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<iostream>
#include<cassert>
#include<string>
namespace material{
struct material_t{
float Z;
float A;
float rho;
float iPot;
float X_0;
};
const material_t Al{ //Konstanten aus Attix
.Z=13.0,
.A=27.0,
.rho=2.69,
.iPot=166.0,
.X_0=24.01
};
const material_t Kapton{ //C22H10N2O5 Konstanten aus Attix
.Z=196.0/39.0,
.A=382.0/39.0,
.rho=1.42,
.iPot=79.6,
.X_0=40.5728 //Attix, Tsai
};
const material_t Air{ //Konstanten aus Leo
.Z=29.0/4.0,
.A=58.0/4.0,
.rho=1.205e-3,
.iPot=85.7,
.X_0=36.66 //Wert aus Tsai
};
const material_t H2O{ //Konstanten aus Leo
.Z=10.0/3.0,
.A=18.0/3.0,
.rho=1.0,
.iPot=75.0,
.X_0=36.08 //Wert aus Tsai
};
const material_t Szintillator{ //Konstanten aus Leo
.Z=7.0/2.0,
.A=13.0/2.0,
.rho=1.06,
.iPot=64.7,
.X_0=43.8
};
const material_t Mylar{ //C5H4O2 Konstanten aus Attix
.Z=50.0/11.0,
.A=96.0/11.0,
.rho=1.40,
.iPot=78.7,
.X_0=39.95
};
const material_t Si{ //Konstanten aus Leo
.Z=14.0,
.A=28.0,
.rho=2.33,
.iPot=173.0,
.X_0=21.82
};
const material_t Au{ //Konstanten aus Leo
.Z=79.0,
.A=197.0,
.rho=19.28,
.iPot=790.0,
.X_0=6.468
};
const material_t Ta{ //Konstanten aus Sternheimer
.Z=73.0,
.A=181.0,
.rho=16.654,
.iPot=718.0,
.X_0=6.8177
};
const material_t Plexiglas{ //Polymethylmethacrulat C5H8O2
//Konstanten aus Sternheimer
.Z=54./15.0,
.A=100./15.0,
.rho=1.19,
.iPot=74.0,
.X_0=40.55
};
const material_t Pb{ //Konstanten aus PDB und Leo
.Z=82.0,
.A=207.2,
.rho=11.35,
.iPot=823.0,
.X_0=6.37
};
const material_t Fe{ //Konstanten aus Leo
.Z=26.0,
.A=55.85,
.rho=7.5,
.iPot=286.0,
.X_0=13.84
};
const material_t LiF{ //Konstanten aus PDB
.Z=(3.+9.)/2.0,
.A=(6.941+18.99840)/2.0,
.rho=2.635,
.iPot=94.0,
.X_0=39.25
};
const material_t TLD{ //LiF Dichte 2.5
//Konstanten aus PDB und Fritsch priv. aus Polen
.Z=(3.+9.)/2.0,
.A=(6.941+18.99840)/2.0,
.rho=2.5, // <-- einziger Unterschied
.iPot=94.0,
.X_0=39.25
};
const material_t PS{ //Konstanten aus Sternheimer
.Z =3.5, //geschätzt
.A =(1.008+12.011)/2.0, // Z/A stimmt
.rho =1.06,
.iPot =68.7,
.X_0 =43.79
};
struct target_t{
material_t mat;
float thicknessInMicrometer;
};
material_t getMaterialFromUser(){
std::cout << "ENERGIEVERLUST NACH BETHE-BLOCH (Formeln aus Leo)" <<std::endl;
std::cout << "-----------------------------------------------------" <<std::endl;
std::cout << "a ___ Energieverlust in Aluminium" <<std::endl;
std::cout << "k ___ Energieverlust in Kapton" <<std::endl;
std::cout << "l ___ Energieverlust in Luft" <<std::endl;
std::cout << "w ___ Energieverlust in Wasser" <<std::endl;
std::cout << "s ___ Energieverlust in Szintillator" <<std::endl;
std::cout << "m ___ Energieverlust in Mylarfolie" <<std::endl;
std::cout << "z ___ Energieverlust in Silizium" <<std::endl;
std::cout << "g ___ Energieverlust in Gold" <<std::endl;
std::cout << "t ___ Energieverlust in Tantal" <<std::endl;
std::cout << "p ___ Energieverlust in Plexiglas" <<std::endl;
std::cout << "b ___ Energieverlust in Blei" <<std::endl;
std::cout << "e ___ Energieverlust in Eisen" <<std::endl;
std::cout << "i ___ Energieverlust in LiF " <<std::endl;
std::cout << "d ___ Energieverlust in TLD (LiF mit Dichte 2.50) " <<std::endl;
std::cout << "r ___ Energieverlust in Polystyrol PS " <<std::endl;
std::cout << "q ___ Ende" <<std::endl;
char selection;
std::cin >> selection;
switch(selection){
default : std::cout << "Eingabe nicht erkannt";
exit(1);
case 'a' : return Al;
case 'k' : return Kapton;
case 'l' : return Air;
case 'w' : return H2O;
case 's' : return Szintillator;
case 'm' : return Mylar;
case 'z' : return Si;
case 'g' : return Au;
case 't' : return Ta;
case 'p' : return Plexiglas;
case 'b' : return Pb;
case 'e' : return Fe;
case 'i' : return LiF;
case 'd' : return TLD;
case 'r' : return PS;
case 'q' : exit(0);
};
}
float getThicknessFromUser(){
std::cout << "Welche Stärke hat das Targetmaterial in µm: " << std::endl;
float thicknessInMicrometer;
std::string inp;
std::cin >> inp;
thicknessInMicrometer = std::stof(inp);
return thicknessInMicrometer;
}
target_t makeTargetFromUser(){
material_t mat = getMaterialFromUser();
float thicknessInMicrometer = getThicknessFromUser();
target_t retval{
.mat = mat,
.thicknessInMicrometer = thicknessInMicrometer
};
return retval;
}
};
namespace dedx{
float getEKinFromUser(){
float Ekin;
std::string inp;
std::cout << "Bitte die Kinetische Energie des Protons eingeben: " << std::endl;
std::cin >> inp;
Ekin = std::stof(inp);
return Ekin;
}
void evaluateInteraction( float EkinProton, material::target_t target){
int i;
float e, impuls,beta,gamma,eta,w_max;
float c1{}, c2{}, c{}, dedx{}, range{};
float EkinProton0{EkinProton};
int zahl{1000};
/* Density-Korrekturen spielen bei den niedrigen Energien noch keine
Rolle => delta = 0 */
float const delta {0.};
const float ProtonMass {938.27231};
const float cnst = 0.1535;
const float thicknessInMillimeter = target.thicknessInMicrometer/1000.0;
for (i=1;i<=zahl;i++) {
e = ProtonMass + EkinProton;
impuls = sqrt(e*e - ProtonMass*ProtonMass);
beta = impuls/e;
int firstBetaLower004 = 0;
if ((beta <= 0.04) && (firstBetaLower004!=0)){
firstBetaLower004 = i;
std::cout << "beta ist ab Iteration " << firstBetaLower004 <<" kleiner als 0.4!"<< std::endl;
std::cout << "Um Fehler zu vermeiden (log negativ etc.) wird beta auf 0.4 festgelegt"
<< std::endl;
beta=0.04;
}
gamma = 1./sqrt(1.-(beta*beta));
eta = beta*gamma;
if (eta >= 0.1)
{
//std::cout << "eta ist größer als 0.1!" << std::endl;
//std::cout << "Shell-Korrektur wird angewendet" << std::endl;
c1 = 0.422377/pow((double)eta,(double)2)
+ 0.0304043/pow((double)eta,(double)4)
- 0.00038106/pow((double)eta,(double)6);
c2 = 3.850190/pow((double)eta,(double)2)
- 0.1667989/pow((double)eta,(double)4)
+ 0.00157955/pow((double)eta,(double)6);
c = c1*1.0e-6*pow((double)target.mat.iPot,(double)2) + c2*1.0e-9
*pow((double)target.mat.iPot,(double)3);
}
else
{
std::cout << "eta ist kleiner als 0.1!" << std::endl;
std::cout << "c wird auf 0 gesetzt um Fehler in der Shell-Korrektur zu vermeiden"
<< std::endl;
c = 0.;
}
dedx = cnst/ (10.0*target.mat.rho*target.mat.Z/target.mat.A*1.0)
/(beta*beta)*(log(2.*0.511*eta*eta*2.*0.511*eta*eta
/(target.mat.iPot*1.e-6*target.mat.iPot*1.e-6))-
2.*beta*beta-delta-2.*c/target.mat.Z);
if (i == 1)
{
std::cout << "Vor der Schicht:" << std::endl;
std::cout << "Kin. Energie (MeV): " << EkinProton << std::endl;
std::cout << "dE/dx (keV/um): " << dedx << std::endl;
std::cout << "Gesamtenergie (MeV): "<< e << std::endl;
std::cout << "Impuls (MeV/c): "<< impuls << std::endl;
std::cout << "Beta: "<< beta << std::endl;
std::cout << "Gamma: "<< gamma << std::endl;
std::cout << "Eta: " << eta << std::endl;
w_max = 2.*0.511*eta*eta/(1.+2.*(0.511/ProtonMass*sqrt(1.+(eta*eta)))
+0.511/ProtonMass);
std::cout << "Max. Energieuebertrag auf ein Elektron (MeV): " << w_max << std::endl;
}
EkinProton = EkinProton - dedx*thicknessInMillimeter;
if (EkinProton <= 0.0) {
range = (float)(i-1)*target.thicknessInMicrometer;
std::cout << "Teilchen bleibt nach ca. "<< range << " um stecken";
exit(0);
}
}
e = ProtonMass + EkinProton;
impuls = sqrt(e*e - ProtonMass*ProtonMass);
beta = impuls/e;
gamma = 1./sqrt(1.-(beta*beta));
eta = beta*gamma;
std::cout << "-----------------------------" << std::endl;
std::cout << "Bei Austritt aus der Schicht:" << std::endl;
std::cout << "Kin. Energie (MeV): " << EkinProton << std::endl;
std::cout << "dE/dx (keV/um): " << dedx << std::endl;
std::cout << "Gesamtenergie (MeV): " << e << std::endl;
std::cout << "Impuls (MeV/c): " << impuls << std::endl;
std::cout << "Beta: "<< beta << std::endl;
std::cout << "Gamma: " << gamma << std::endl;
std::cout << "Eta: " << eta << std::endl;
w_max = 2.*0.511*eta*eta/(1.+2.*(0.511/ProtonMass*sqrt(1.+(eta*eta)))
+0.511/ProtonMass);
std::cout << "Max. Energieuebertrag auf ein Elektron (MeV): " << w_max << std::endl;
std::cout << "Energieverlust in der Schicht (MeV): " << (EkinProton0 - EkinProton) << std::endl;
if (beta <= 0.04)
std::cout << "ACHTUNG: beta kleiner 0.04, Ergebnis nach Bethe-Bloch-Formel ist falsch!" << std::endl;
return;
}
};
dedx.cpp 0 → 100644
+ 13
− 0
View file @ 00922e02
#include "bethebloch.hpp"
#include "iostream"
int main(){
material::target_t target = material::makeTargetFromUser();
float EkinProton = dedx::getEKinFromUser();
dedx::evaluateInteraction(EkinProton, target);
}
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment