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 12a89c9..54a4414 100644 --- a/input.txt +++ b/input.txt @@ -1,26 +1,32 @@ ----------Data Information---------- -OutputFile: /data1/gwm17/test_dead.root +OutputFile: /media/gordon/b6414c35-ec1f-4fc1-83bc-a6b68ca4325a/gwm17/test_newkine.root SaveTree: yes SavePlots: yes ----------Reaction Information---------- ReactionType: 2 Z A (order is target, projectile, ejectile, break1, break3) -6 12 +5 10 2 3 +2 4 1 2 -1 1 ----------Target Information---------- Name: test_targ -Layers: 1 +Layers: 2 ~Layer1 -Thickness(ug/cm^2): 40 +Thickness(ug/cm^2): 10 Z A Stoich 6 12 1 0 ~ +~Layer2 +Thickness(ug/cm^2): 80 +Z A Stoich +5 10 1 +0 +~ ----------Sampling Information---------- NumberOfSamples: 1000000 BeamMeanEnergy(MeV): 24 BeamEnergySigma(MeV): 0.001 -EjectileThetaMin(deg): 3.0 EjectileThetaMax(deg): 3.0 -ResidualExMean(MeV): 2.364 ResidualExSigma(MeV): 0.0317 +EjectileThetaMin(deg): 20.0 EjectileThetaMax(deg): 20.0 +ResidualExMean(MeV): 16.8 ResidualExSigma(MeV): 0.038 -------------------------------------- diff --git a/input.txt.orig b/input.txt.orig new file mode 100644 index 0000000..72f9646 --- /dev/null +++ b/input.txt.orig @@ -0,0 +1,46 @@ +----------Data Information---------- +OutputFile: /data1/gwm17/test_dead.root +SaveTree: yes +SavePlots: yes +----------Reaction Information---------- +ReactionType: 2 +Z A (order is target, projectile, ejectile, break1, break3) +<<<<<<< HEAD +6 12 +2 3 +1 2 +1 1 +======= +5 10 +2 3 +2 4 +1 2 +>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7 +----------Target Information---------- +Name: test_targ +Layers: 1 +~Layer1 +<<<<<<< HEAD +Thickness(ug/cm^2): 40 +======= +Thickness(ug/cm^2): 10 +>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7 +Z A Stoich +6 12 1 +0 +~ +<<<<<<< HEAD +======= +~Layer2 +Thickness(ug/cm^2): 80 +Z A Stoich +5 10 1 +0 +~ +>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7 +----------Sampling Information---------- +NumberOfSamples: 1000000 +BeamMeanEnergy(MeV): 24 BeamEnergySigma(MeV): 0.001 +EjectileThetaMin(deg): 3.0 EjectileThetaMax(deg): 3.0 +ResidualExMean(MeV): 2.364 ResidualExSigma(MeV): 0.0317 +-------------------------------------- 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(); diff --git a/src/Target.cpp b/src/Target.cpp index 99dc473..24b5b35 100644 --- a/src/Target.cpp +++ b/src/Target.cpp @@ -40,28 +40,28 @@ 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 if(theta > PI/2.) theta = PI - 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 if(theta > PI/2.) theta = PI - 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 if(theta > PI/2.) theta = PI - 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 if(theta > PI/2.) theta = PI - theta; + else if (theta > PI/2.) theta = PI - theta; return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta)))); } diff --git a/src/Target.cpp.orig b/src/Target.cpp.orig new file mode 100644 index 0000000..8e47272 --- /dev/null +++ b/src/Target.cpp.orig @@ -0,0 +1,104 @@ +/* + +Target.cpp +A basic target unit for use in the SPANCRedux environment. A target +is defined as a single compound with elements Z,A of a given stoichiometry +Holds an energy loss class + +Based on code by D.W. Visser written at Yale for the original SPANC + +Written by G.W. McCann Aug. 2020 + +*/ +#include "Target.h" + +/*Targets must be of known thickness*/ +Target::Target(double thick) { + thickness = thick; +} + +Target::~Target() { +} + +/*Set target elements of given Z, A, S*/ +void Target::SetElements(std::vector& z, std::vector& a, std::vector& stoich) { + Z = z; + A = a; + Stoich = stoich; + + eloss.SetTargetComponents(Z, A, Stoich); +} + +/*Element verification*/ +bool Target::ContainsElement(int z, int a) { + for(unsigned int i=0; i PI/2.) theta = PI - theta; +======= + else if (theta > PI/2.) theta = PI - theta; +>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7 + 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; +<<<<<<< HEAD + else if(theta > PI/2.) theta = PI - theta; +======= + else if (theta > PI/2.) theta = PI - theta; +>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7 + 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; +<<<<<<< HEAD + else if(theta > PI/2.) theta = PI - theta; +======= + else if (theta > PI/2.) theta = PI - theta; +>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7 + 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; +<<<<<<< HEAD + else if(theta > PI/2.) theta = PI - theta; +======= + else if (theta > PI/2.) theta = PI - theta; +>>>>>>> 4ccaabb534f35f2cf36d375b74ac4ebf99fe7bb7 + return eloss.GetReverseEnergyLoss(zp, ap, finalEnergy, thickness/(2.0*fabs(cos(theta)))); +} + +/*Getter functions*/ + +double& Target::GetThickness() { + return thickness; +} + +int Target::GetNumberOfElements() { + return Z.size(); +} + +int Target::GetElementZ(int index) { + return Z[index]; +} + +int Target::GetElementA(int index) { + return A[index]; +} + +int Target::GetElementStoich(int index) { + return Stoich[index]; +}