From dc7f3923edb3a5f26e4396032bfd9394c148fcee Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Sat, 13 Feb 2021 17:35:29 -0500 Subject: [PATCH 1/2] Fixed bug in target thickness calc, where angle through target was incorrectly used for energy loss --- src/Target.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/Target.cpp b/src/Target.cpp index fffd480..b3757df 100644 --- a/src/Target.cpp +++ b/src/Target.cpp @@ -40,25 +40,29 @@ bool Target::ContainsElement(int z, int a) { /*Calculates energy loss for travelling all the way through the target*/ double Target::getEnergyLossTotal(int zp, int ap, double startEnergy, double theta) { if(theta == PI/2.) return startEnergy; - else return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta))); + else if (theta > PI/2.) theta = PI - theta; + return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/fabs(cos(theta))); } /*Calculates energy loss for travelling halfway through the target*/ double Target::getEnergyLossHalf(int zp, int ap, double startEnergy, double theta) { if(theta == PI/2.) return startEnergy; - else return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/(2.0*fabs(cos(theta)))); + else if (theta > PI/2.) theta = PI - theta; + return eloss.GetEnergyLoss(zp, ap, startEnergy, thickness/(2.0*fabs(cos(theta)))); } /*Calculates reverse energy loss for travelling all the way through the target*/ double Target::getReverseEnergyLossTotal(int zp, int ap, double finalEnergy, double theta) { if(theta == PI/2.) return finalEnergy; - else return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta))); + else if (theta > PI/2.) theta = PI - theta; + return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/fabs(cos(theta))); } /*Calculates reverse energy loss for travelling half way through the target*/ double Target::getReverseEnergyLossHalf(int zp, int ap, double finalEnergy, double theta) { if(theta == PI/2.) return finalEnergy; - else return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta)))); + else if (theta > PI/2.) theta = PI - theta; + return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta)))); } /*Getter functions*/ From 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7 Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Sat, 13 Feb 2021 17:38:38 -0500 Subject: [PATCH 2/2] Some commenting added --- include/EnergyLoss.h | 1 + include/Nucleus.h | 11 ++++++++++- include/Reaction.h | 10 ++++++++++ include/ReactionSystem.h | 14 ++++++++++++++ input.txt | 15 ++++++++------- src/EnergyLoss.cpp | 8 +++++++- src/Nucleus.cpp | 10 +++++++++- src/Reaction.cpp | 10 ++++++++++ 8 files changed, 69 insertions(+), 10 deletions(-) diff --git a/include/EnergyLoss.h b/include/EnergyLoss.h index 4ae90ae..fc0283f 100644 --- a/include/EnergyLoss.h +++ b/include/EnergyLoss.h @@ -46,6 +46,7 @@ class EnergyLoss { //constants for calculations static constexpr double MAX_FRACTIONAL_STEP = 0.001; + static constexpr double MAX_DEPTH = 50; static constexpr double MAX_H_E_PER_U = 100000.0; static constexpr double AVOGADRO = 0.60221367; //N_A times 10^(-24) for converting static constexpr double MEV2U = 1.0/931.4940954; diff --git a/include/Nucleus.h b/include/Nucleus.h index 548c200..2fd0c64 100644 --- a/include/Nucleus.h +++ b/include/Nucleus.h @@ -1,3 +1,11 @@ +/* + Nucleus.h + Nucleus is a derived class of Vec4. A nucleus is the kinematics is essentially a 4 vector with the + additional properties of the number of total nucleons (A), the number of protons (Z), a ground state mass, + an exctitation energy, and an isotopic symbol. + + --GWM Jan 2021 +*/ #ifndef NUCLEUS_H #define NUCLEUS_H @@ -25,7 +33,8 @@ public: return *this; }; - inline Nucleus operator+(const Nucleus& daughter) { + //Conservation of nucleons and momentum + inline Nucleus operator+(const Nucleus& daughter) { return Nucleus(GetZ()+daughter.GetZ(), GetA()+daughter.GetA(), GetPx()+daughter.GetPx(), GetPy()+daughter.GetPy(), GetPz()+daughter.GetPz(), GetE()+daughter.GetE()); }; inline Nucleus operator-(const Nucleus& daughter) { diff --git a/include/Reaction.h b/include/Reaction.h index 7b72ad1..d8b1b40 100644 --- a/include/Reaction.h +++ b/include/Reaction.h @@ -1,3 +1,11 @@ +/* + Reaction.h + Reaction is a class which implements either a decay or scattering reaction. As such it requires either + 3 (decay) or 4 (scattering) nuclei to perform any calcualtions. I also links together the target, which provides + energy loss calculations, with the kinematics. Note that Reaction does not own the LayeredTarget. + + --GWM Jan. 2021 +*/ #ifndef REACTION_H #define REACTION_H @@ -15,6 +23,8 @@ public: void SetNuclei(int zt, int at, int zp, int ap, int ze, int ae); void SetNuclei(const Nucleus* nucs); void SetBeamKE(double bke); + + /*Setters and getters*/ inline void SetLayeredTarget(LayeredTarget* targ) { target = targ; }; inline void SetPolarRxnAngle(double theta) { m_theta = theta; }; inline void SetAzimRxnAngle(double phi) { m_phi = phi; }; diff --git a/include/ReactionSystem.h b/include/ReactionSystem.h index 6f86739..63fef79 100644 --- a/include/ReactionSystem.h +++ b/include/ReactionSystem.h @@ -1,3 +1,11 @@ +/* + ReactionSystem.h + ReactionSystem is the base class from which all kinematics calculations should inherit. It contains all members and functions + to perform a single step reaction/decay, which is the basic building block of all subsequent types of reactions in MASK. More + complicated systems (see TwoStepSystem and ThreeStepSystem) add further data members and override the virtual functions. + + --GWM Jan. 2021 +*/ #ifndef REACTIONSYSTEM_H #define REACTIONSYSTEM_H @@ -12,12 +20,16 @@ public: ReactionSystem(); ReactionSystem(std::vector& z, std::vector& a); virtual ~ReactionSystem(); + virtual bool SetNuclei(std::vector& z, std::vector& a); void AddTargetLayer(std::vector& zt, std::vector& at, std::vector& stoich, double thickness); + + /*Set sampling parameters*/ inline void SetRandomGenerator(TRandom3* gen) { generator = gen; gen_set_flag = true; }; inline void SetBeamDistro(double mean, double sigma) { m_beamDist = std::make_pair(mean, sigma); }; inline void SetTheta1Range(double min, double max) { m_theta1Range = std::make_pair(min*deg2rad, max*deg2rad); }; inline void SetExcitationDistro(double mean, double sigma) { m_exDist = std::make_pair(mean, sigma); }; + virtual void RunSystem(); inline const Nucleus& GetTarget() const { return step1.GetTarget(); }; @@ -32,6 +44,8 @@ protected: Reaction step1; LayeredTarget target; + + //Sampling information std::pair m_beamDist, m_theta1Range, m_exDist; TRandom3* generator; //not owned by ReactionSystem diff --git a/input.txt b/input.txt index 0dfa5e0..54a4414 100644 --- a/input.txt +++ b/input.txt @@ -3,24 +3,25 @@ OutputFile: /media/gordon/b6414c35-ec1f-4fc1-83bc-a6b68ca4325a/gwm17/test_newkin SaveTree: yes SavePlots: yes ----------Reaction Information---------- -ReactionType: 0 +ReactionType: 2 Z A (order is target, projectile, ejectile, break1, break3) -5 9 -0 0 -1 1 +5 10 +2 3 +2 4 +1 2 ----------Target Information---------- Name: test_targ Layers: 2 ~Layer1 -Thickness(ug/cm^2): 0 +Thickness(ug/cm^2): 10 Z A Stoich 6 12 1 0 ~ ~Layer2 -Thickness(ug/cm^2): 0 +Thickness(ug/cm^2): 80 Z A Stoich -5 9 1 +5 10 1 0 ~ ----------Sampling Information---------- diff --git a/src/EnergyLoss.cpp b/src/EnergyLoss.cpp index a5a9e19..4e970c1 100644 --- a/src/EnergyLoss.cpp +++ b/src/EnergyLoss.cpp @@ -55,6 +55,8 @@ double EnergyLoss::GetEnergyLoss(int zp, int ap, double e_initial, double thickn double e_step = GetTotalStoppingPower(e_final)*x_step/1000.0; //initial step in e double e_threshold = 0.05*e_initial; + int depth=0; + if(thickness == 0.0) return 0; else if(e_initial == 0.0) return 0; @@ -62,16 +64,20 @@ double EnergyLoss::GetEnergyLoss(int zp, int ap, double e_initial, double thickn bool go = true; while(go) { //If intial guess of step size is too large, shrink until in range - if(e_step/e_final > MAX_FRACTIONAL_STEP /*&& e_step >= E_PRECISION_LIMIT*/) { + if(e_step/e_final > MAX_FRACTIONAL_STEP && depth < MAX_DEPTH) { + depth++; x_step *= 0.5; e_step = GetTotalStoppingPower(e_final)*x_step/1000.0; } else if((x_step + x_traversed) >= thickness) { //last chunk go = false; x_step = thickness - x_traversed; //get valid portion of last chunk e_final -= GetTotalStoppingPower(e_final)*x_step/1000.0; + if(depth > 20)std::cout<<"depth: "<GetProjectileEnergyLoss(reactants[1].GetZ(), reactants[1].GetA(), bke, rxnLayer, 0); }; +//Methods given by Iliadis in Nuclear Physics of Stars, Appendix C void Reaction::CalculateReaction() { //Target assumed at rest, with 0 excitation energy reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass()); @@ -112,6 +121,7 @@ void Reaction::CalculateReaction() { } +//Calculate in CM, where decay is isotropic void Reaction::CalculateDecay() { double Q = reactants[0].GetInvMass() - reactants[2].GetGroundStateMass() - reactants[3].GetGroundStateMass();