diff --git a/src/Detectors/SabreArray.cpp b/src/Detectors/SabreArray.cpp index 49c9a01..cef9093 100644 --- a/src/Detectors/SabreArray.cpp +++ b/src/Detectors/SabreArray.cpp @@ -20,6 +20,12 @@ SabreArray::SabreArray() : m_degradedDetectors[2] = false; m_degradedDetectors[3] = false; m_degradedDetectors[4] = true; + //No degraded detectors + // m_degradedDetectors[0] = false; + // m_degradedDetectors[1] = false; + // m_degradedDetectors[2] = false; + // m_degradedDetectors[3] = false; + // m_degradedDetectors[4] = false; //Choose who to look at right now. Usually switch on or off degraded/non-degraded. m_activeDetectors[0] = false; @@ -142,4 +148,4 @@ DetectorResult SabreArray::IsDetected(const Mask::Nucleus& nucleus) observation.detectFlag = false; return observation; -} \ No newline at end of file +} diff --git a/src/Mask/AngularDistribution.cpp b/src/Mask/AngularDistribution.cpp index 8f07b2b..b21ddf0 100644 --- a/src/Mask/AngularDistribution.cpp +++ b/src/Mask/AngularDistribution.cpp @@ -79,7 +79,7 @@ namespace Mask { //Renormalize distribution such that total prob is 1.0. //Test branching ratio to see if we "make" a decay particle, //then use re-normalized distribution to pick an angle. - if(m_constants[0] < 0.5) + if(m_constants[0] != 0.5) { double norm = 0.5/m_constants[0]; for(auto& value : m_constants) @@ -115,4 +115,14 @@ namespace Mask { return costheta; } + double AngularDistribution::GetProbability(double cosTheta) + { + double prob = 0.0; + for (std::size_t i=0; i& params) : + ReactionSystem() + { + m_nuclei.resize(8); + Init(params); + } + + CoupledThreeStepSystem::~CoupledThreeStepSystem() {} + + void CoupledThreeStepSystem::Init(const std::vector& params) + { + if(params.size() != 3 || params[0].rxnType != RxnType::Reaction || params[1].rxnType != RxnType::Decay || + params[2].rxnType != RxnType::Decay || params[0].Z.size() != 3 || params[0].A.size() != 3 || params[1].Z.size() != 2 || + params[1].A.size() != 2 || params[2].Z.size() != 2 || params[2].A.size() != 2) + { + m_isValid = false; + std::cerr << "Invalid parameters at ThreeStepSystem::Init(), does not match ThreeStep signature!" << std::endl; + return; + } + + const StepParameters& step1Params = params[0]; + const StepParameters& step2Params = params[1]; + const StepParameters& step3Params = params[2]; + + //Setup nuclei + int zr = step1Params.Z[0] + step1Params.Z[1] - step1Params.Z[2]; + int ar = step1Params.A[0] + step1Params.A[1] - step1Params.A[2]; + if(zr != step2Params.Z[0] || ar != step2Params.A[0]) + { + m_isValid = false; + std::cerr << "Invalid parameters at ThreeStepSystem::Init(), step one and step two are not sequential! Step one recoil (Z,A): (" + << zr << "," << ar << ") Step two target (Z,A): (" << step2Params.Z[0] << "," << step2Params.A[0] << ")" <" + << m_nuclei[4].isotopicSymbol << "+" + << m_nuclei[5].isotopicSymbol << "->" + << m_nuclei[6].isotopicSymbol << "+" + << m_nuclei[7].isotopicSymbol; + m_sysEquation = stream.str(); + } + + CoupledThreeStepParameters CoupledThreeStepSystem::SampleParameters() + { + CoupledThreeStepParameters params; + std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); + params.beamEnergy = m_beamDistributions[0](gen); + params.rxnTheta = std::acos((m_thetaRanges[0])(gen)); + params.rxnPhi = m_phiRanges[0](gen); + params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta(); + params.cosRelativeAngle = m_decayAngularDistributions[1].GetRandomCosTheta(); + params.decay1Theta = std::acos(params.cosdecay1Theta); + params.decay1Phi = m_phiRanges[1](gen); + params.residEx = m_exDistributions[0](gen); + params.decay1Ex = m_exDistributions[1](gen); + params.decay2Ex = m_exDistributions[2](gen); + params.rxnDepth = (m_rxnDepthDist(gen)); + return params; + } + + //Called after running step2 + void CoupledThreeStepSystem::SampleCoupling(CoupledThreeStepParameters& params) + { + ROOT::Math::Boost boost(m_nuclei[5].vec4.BoostToCM()); + ROOT::Math::PxPyPzEVector a1Vec = (boost * m_nuclei[4].vec4); + ROOT::Math::XYZVector a1PVec(a1Vec.Px(), a1Vec.Py(), a1Vec.Pz()); + ROOT::Math::PxPyPzEVector a2Vec; + ROOT::Math::XYZVector a2PVec; + ROOT::Math::RotationZ zRot; + ROOT::Math::RotationY yRot; + a1PVec *= 1.0/a1Vec.P(); + zRot.SetAngle(a1PVec.Phi()); + yRot.SetAngle(a1PVec.Theta()); + double relAngle = std::acos(params.cosRelativeAngle); + ROOT::Math::XYZVector a1PVecAligned = yRot.Inverse() * (zRot.Inverse() * a1PVec); + a2PVec.SetXYZ( + std::sin(relAngle), + 0.0, + std::cos(relAngle) + ); + a2PVec = zRot * (yRot * a2PVec); + params.decay2Theta = a2PVec.Theta(); + params.decay2Phi = a2PVec.Phi(); + } + + void CoupledThreeStepSystem::RunSystem() + { + CoupledThreeStepParameters params; + do + { + params = SampleParameters(); + } + while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) + && m_step2.CheckDecayThreshold(params.residEx, params.decay1Ex) + && m_step3.CheckDecayThreshold(params.decay1Ex, params.decay2Ex))); + + m_step1.SetReactionDepth(params.rxnDepth); + m_step1.SetBeamKE(params.beamEnergy); + m_step1.SetPolarRxnAngle(params.rxnTheta); + m_step1.SetAzimRxnAngle(params.rxnPhi); + m_step1.SetExcitation(params.residEx); + + m_step2.SetReactionDepth(params.rxnDepth); + m_step2.SetPolarRxnAngle(params.decay1Theta); + m_step2.SetAzimRxnAngle(params.decay1Phi); + m_step2.SetExcitation(params.decay1Ex); + + m_step3.SetReactionDepth(params.rxnDepth); + + m_step3.SetExcitation(params.decay2Ex); + + m_step1.Calculate(); + + if(params.cosdecay1Theta == -10) + { + m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass); + m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass); + m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass); + m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass); + return; + } + m_step2.Calculate(); + + m_step3.SetResidualEnergyLoss(true); + + SampleCoupling(params); //must be called here as frame needs to be set + + m_step3.SetPolarRxnAngle(params.decay2Theta); + m_step3.SetAzimRxnAngle(params.decay2Phi); + m_step3.Calculate(); + } +} \ No newline at end of file diff --git a/src/Mask/CoupledThreeStepSystem.h b/src/Mask/CoupledThreeStepSystem.h new file mode 100644 index 0000000..fed700e --- /dev/null +++ b/src/Mask/CoupledThreeStepSystem.h @@ -0,0 +1,46 @@ +#ifndef COUPLEDTHREESTEPSYSTEM_H +#define COUPLEDTHREESTEPSYSTEM_H + +#include "ReactionSystem.h" +#include "AngularDistribution.h" + +namespace Mask { + + struct CoupledThreeStepParameters + { + double beamEnergy = 0.; + double rxnTheta = 0.; + double rxnPhi = 0.; + double cosdecay1Theta = 0.; + double cosRelativeAngle = 0.; + double decay1Theta = 0.; + double decay1Phi = 0.; + double decay2Theta = 0.; + double decay2Phi = 0.; + double residEx = 0.; + double decay1Ex = 0.; + double decay2Ex = 0.; + double rxnDepth = 0.; + }; + + class CoupledThreeStepSystem : public ReactionSystem + { + public: + CoupledThreeStepSystem(const std::vector& params); + ~CoupledThreeStepSystem(); + + virtual void SetLayeredTarget(const LayeredTarget& target) override; + virtual void RunSystem() override; + + protected: + void Init(const std::vector& params); + void SetSystemEquation() override; + CoupledThreeStepParameters SampleParameters(); + void SampleCoupling(CoupledThreeStepParameters& params); + + Reaction m_step1, m_step2, m_step3; + }; + +} + +#endif \ No newline at end of file diff --git a/src/Mask/Reaction.cpp b/src/Mask/Reaction.cpp index 2a1a7f8..1a1b871 100644 --- a/src/Mask/Reaction.cpp +++ b/src/Mask/Reaction.cpp @@ -256,6 +256,27 @@ namespace Mask { } } + bool Reaction::CheckReactionThreshold(double beamEnergy, double residEx) + { + double Q = m_target->groundStateMass + m_projectile->groundStateMass - + (m_ejectile->groundStateMass + m_residual->groundStateMass + residEx); + double Ethresh = -Q*(m_ejectile->groundStateMass+m_residual->groundStateMass) / + (m_ejectile->groundStateMass + m_residual->groundStateMass - m_projectile->groundStateMass); + if (beamEnergy < Ethresh) + return false; + else + return true; + } + + bool Reaction::CheckDecayThreshold(double targetEx, double residEx) + { + double Q = m_target->groundStateMass + targetEx - + (m_ejectile->groundStateMass + m_residual->groundStateMass + residEx); + if ( Q < 0.0) + return false; + else + return true; + } } diff --git a/src/Mask/Reaction.h b/src/Mask/Reaction.h index 12b9e07..f6ea1e3 100644 --- a/src/Mask/Reaction.h +++ b/src/Mask/Reaction.h @@ -42,6 +42,9 @@ namespace Mask { void SetRxnLayer(std::size_t layer) { m_rxnLayer = layer; }; void SetResidualEnergyLoss(bool isEloss) { m_isResidEloss = isEloss; }; + bool CheckReactionThreshold(double beamEnergy, double residualEx); + bool CheckDecayThreshold(double targetEx, double residualEx); + bool IsDecay() const { return m_isDecay; }; std::size_t GetRxnLayer() const { return m_rxnLayer; }; diff --git a/src/Mask/ThreeStepSystem.cpp b/src/Mask/ThreeStepSystem.cpp index 65e1206..4d43acb 100644 --- a/src/Mask/ThreeStepSystem.cpp +++ b/src/Mask/ThreeStepSystem.cpp @@ -2,6 +2,11 @@ #include "RandomGenerator.h" #include "KinematicsExceptions.h" +#include "Math/Boost.h" +#include "Math/Vector3D.h" +#include "Math/RotationZ.h" +#include "Math/RotationY.h" + namespace Mask { ThreeStepSystem::ThreeStepSystem(const std::vector& params) : @@ -114,44 +119,58 @@ namespace Mask { << m_nuclei[7].isotopicSymbol; m_sysEquation = stream.str(); } + + ThreeStepParameters ThreeStepSystem::SampleParameters() + { + ThreeStepParameters params; + std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); + params.beamEnergy = m_beamDistributions[0](gen); + params.rxnTheta = std::acos((m_thetaRanges[0])(gen)); + params.rxnPhi = m_phiRanges[0](gen); + params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta(); + params.decay1Theta = std::acos(params.cosdecay1Theta); + params.decay1Phi = m_phiRanges[1](gen); + params.cosdecay2Theta = m_decayAngularDistributions[1].GetRandomCosTheta(); + params.decay2Theta = std::acos(params.cosdecay2Theta); + params.decay2Phi = m_phiRanges[1](gen); + params.residEx = m_exDistributions[0](gen); + params.decay1Ex = m_exDistributions[1](gen); + params.decay2Ex = m_exDistributions[2](gen); + params.rxnDepth = (m_rxnDepthDist(gen)); + return params; + } void ThreeStepSystem::RunSystem() { - //Sample parameters - std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); - double bke = m_beamDistributions[0](gen); - double rxnTheta = std::acos((m_thetaRanges[0])(gen)); - double rxnPhi = m_phiRanges[0](gen); - double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta(); - double decay1Theta = std::acos(decay1costheta); - double decay1Phi = m_phiRanges[1](gen); - double decay2costheta = m_decayAngularDistributions[1].GetRandomCosTheta(); - double decay2Theta = std::acos(decay2costheta); - double decay2Phi = m_phiRanges[2](gen); - double residEx = m_exDistributions[0](gen); - double decay1Ex = m_exDistributions[1](gen); - double decay2Ex = m_exDistributions[2](gen); - double rxnDepth = (m_rxnDepthDist(gen)); + ThreeStepParameters params; + do + { + params = SampleParameters(); + } + while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) + && m_step2.CheckDecayThreshold(params.residEx, params.decay1Ex) + && m_step3.CheckDecayThreshold(params.decay1Ex, params.decay2Ex))); - m_step1.SetReactionDepth(rxnDepth); - m_step1.SetBeamKE(bke); - m_step1.SetPolarRxnAngle(rxnTheta); - m_step1.SetAzimRxnAngle(rxnPhi); - m_step1.SetExcitation(residEx); + m_step1.SetReactionDepth(params.rxnDepth); + m_step1.SetBeamKE(params.beamEnergy); + m_step1.SetPolarRxnAngle(params.rxnTheta); + m_step1.SetAzimRxnAngle(params.rxnPhi); + m_step1.SetExcitation(params.residEx); - m_step2.SetReactionDepth(rxnDepth); - m_step2.SetPolarRxnAngle(decay1Theta); - m_step2.SetAzimRxnAngle(decay1Phi); - m_step2.SetExcitation(decay1Ex); + m_step2.SetReactionDepth(params.rxnDepth); + m_step2.SetPolarRxnAngle(params.decay1Theta); + m_step2.SetAzimRxnAngle(params.decay1Phi); + m_step2.SetExcitation(params.decay1Ex); - m_step3.SetReactionDepth(rxnDepth); - m_step3.SetPolarRxnAngle(decay2Theta); - m_step3.SetAzimRxnAngle(decay2Phi); - m_step3.SetExcitation(decay2Ex); + m_step3.SetReactionDepth(params.rxnDepth); + m_step3.SetPolarRxnAngle(params.decay2Theta); + m_step3.SetAzimRxnAngle(params.decay2Phi); + m_step3.SetResidualEnergyLoss(true); + m_step3.SetExcitation(params.decay2Ex); m_step1.Calculate(); - if(decay1costheta == -10) + if(params.cosdecay1Theta == -10) { m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass); m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass); @@ -161,15 +180,13 @@ namespace Mask { } m_step2.Calculate(); - if(decay2costheta == -10) + if(params.cosdecay2Theta == -10) { m_nuclei[6].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[6].groundStateMass); m_nuclei[7].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[7].groundStateMass); return; } - m_step3.SetResidualEnergyLoss(true); m_step3.Calculate(); - } } \ No newline at end of file diff --git a/src/Mask/ThreeStepSystem.h b/src/Mask/ThreeStepSystem.h index 7b933f4..5238d88 100644 --- a/src/Mask/ThreeStepSystem.h +++ b/src/Mask/ThreeStepSystem.h @@ -6,6 +6,23 @@ namespace Mask { + struct ThreeStepParameters + { + double beamEnergy = 0.; + double rxnTheta = 0.; + double rxnPhi = 0.; + double cosdecay1Theta = 0.; + double cosdecay2Theta = 0.; + double decay1Theta = 0.; + double decay1Phi = 0.; + double decay2Theta = 0.; + double decay2Phi = 0.; + double residEx = 0.; + double decay1Ex = 0.; + double decay2Ex = 0.; + double rxnDepth = 0.; + }; + class ThreeStepSystem : public ReactionSystem { public: @@ -18,6 +35,7 @@ namespace Mask { protected: void Init(const std::vector& params); void SetSystemEquation() override; + ThreeStepParameters SampleParameters(); Reaction m_step1, m_step2, m_step3; diff --git a/src/Mask/TwoStepSystem.cpp b/src/Mask/TwoStepSystem.cpp index 691367f..2b0fcc3 100644 --- a/src/Mask/TwoStepSystem.cpp +++ b/src/Mask/TwoStepSystem.cpp @@ -92,35 +92,59 @@ namespace Mask { << m_nuclei[5].isotopicSymbol; m_sysEquation = stream.str(); } + + TwoStepParameters TwoStepSystem::SampleParameters() + { + TwoStepParameters params; + std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); + params.beamEnergy = (m_beamDistributions[0])(gen); + params.rxnTheta = std::acos((m_thetaRanges[0])(gen)); + params.rxnPhi = (m_phiRanges[0])(gen); + params.cosdecay1Theta = m_decayAngularDistributions[0].GetRandomCosTheta(); + params.decay1Theta = std::acos(params.cosdecay1Theta); + params.decay1Phi = m_phiRanges[1](gen); + params.residEx = (m_exDistributions[0])(gen); + params.decay2Ex = m_exDistributions[1](gen); + params.rxnDepth = (m_rxnDepthDist(gen)); + return params; + } void TwoStepSystem::RunSystem() { //Sample parameters - std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); - double bke = (m_beamDistributions[0])(gen); - double rxnTheta = std::acos((m_thetaRanges[0])(gen)); - double rxnPhi = (m_phiRanges[0])(gen); - double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta(); - double decay1Theta = std::acos(decay1costheta); - double decay1Phi = m_phiRanges[1](gen); - double residEx = (m_exDistributions[0])(gen); - double decay2Ex = m_exDistributions[1](gen); - double rxnDepth = (m_rxnDepthDist(gen)); + // std::mt19937& gen = RandomGenerator::GetInstance().GetGenerator(); + // double bke = (m_beamDistributions[0])(gen); + // double rxnTheta = std::acos((m_thetaRanges[0])(gen)); + // double rxnPhi = (m_phiRanges[0])(gen); + // double decay1costheta = m_decayAngularDistributions[0].GetRandomCosTheta(); + // double decay1Theta = std::acos(decay1costheta); + // double decay1Phi = m_phiRanges[1](gen); + // double residEx = (m_exDistributions[0])(gen); + // double decay2Ex = m_exDistributions[1](gen); + // double rxnDepth = (m_rxnDepthDist(gen)); + + TwoStepParameters params; + do + { + params = SampleParameters(); + } + while(!(m_step1.CheckReactionThreshold(params.beamEnergy, params.residEx) + && m_step2.CheckDecayThreshold(params.residEx, params.decay2Ex))); - m_step1.SetReactionDepth(rxnDepth); - m_step1.SetBeamKE(bke); - m_step1.SetPolarRxnAngle(rxnTheta); - m_step1.SetAzimRxnAngle(rxnPhi); - m_step1.SetExcitation(residEx); + m_step1.SetReactionDepth(params.rxnDepth); + m_step1.SetBeamKE(params.beamEnergy); + m_step1.SetPolarRxnAngle(params.rxnTheta); + m_step1.SetAzimRxnAngle(params.rxnPhi); + m_step1.SetExcitation(params.residEx); - m_step2.SetReactionDepth(rxnDepth); - m_step2.SetPolarRxnAngle(decay1Theta); - m_step2.SetAzimRxnAngle(decay1Phi); - m_step2.SetExcitation(decay2Ex); + m_step2.SetReactionDepth(params.rxnDepth); + m_step2.SetPolarRxnAngle(params.decay1Theta); + m_step2.SetAzimRxnAngle(params.decay1Phi); + m_step2.SetExcitation(params.decay2Ex); m_step1.Calculate(); - if(decay1costheta == -10) + if(params.cosdecay1Theta == -10) { m_nuclei[4].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[4].groundStateMass); m_nuclei[5].vec4.SetPxPyPzE(0., 0., 0., m_nuclei[5].groundStateMass); diff --git a/src/Mask/TwoStepSystem.h b/src/Mask/TwoStepSystem.h index 70a1773..16d29df 100644 --- a/src/Mask/TwoStepSystem.h +++ b/src/Mask/TwoStepSystem.h @@ -6,6 +6,19 @@ namespace Mask { + struct TwoStepParameters + { + double beamEnergy = 0.; + double rxnTheta = 0.; + double rxnPhi = 0.; + double cosdecay1Theta = 0.; + double decay1Theta = 0.; + double decay1Phi = 0.; + double residEx = 0.; + double decay2Ex = 0.; + double rxnDepth = 0.; + }; + class TwoStepSystem : public ReactionSystem { public: @@ -18,6 +31,7 @@ namespace Mask { private: void Init(const std::vector& params); void SetSystemEquation() override; + TwoStepParameters SampleParameters(); Reaction m_step1, m_step2; }; diff --git a/src/Plotters/RootPlotter.cpp b/src/Plotters/RootPlotter.cpp index 1b56ab3..a56b8ae 100644 --- a/src/Plotters/RootPlotter.cpp +++ b/src/Plotters/RootPlotter.cpp @@ -1,7 +1,8 @@ #include "RootPlotter.h" #include #include -#include +#include +#include "Math/Boost.h" #include @@ -10,7 +11,8 @@ static double FullPhi(double phi) return phi < 0.0 ? 2.0*M_PI + phi : phi; } -RootPlotter::RootPlotter() +RootPlotter::RootPlotter() : + m_target({5},{9},{1}, 74.0) { TH1::AddDirectory(kFALSE); } @@ -42,10 +44,16 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam std::cout<<"\rPercent of data processed: "<GetEntry(i); - for(Mask::Nucleus& nuc : *(dataHandle)) + // for(Mask::Nucleus& nuc : *(dataHandle)) + // { + // FillData(nuc); + // } + for(int i=0; isize(); i++) { - FillData(nuc); + FillData(dataHandle->at(i), i); } + //Don't leave this in! + //Correlations(*dataHandle); } std::cout<Close(); @@ -57,34 +65,35 @@ void RootPlotter::Run(const std::string& inputname, const std::string& outputnam output->Close(); } -void RootPlotter::FillData(const Mask::Nucleus& nuc) +void RootPlotter::FillData(const Mask::Nucleus& nuc, int i) { - std::string mod = "detected"; + std::string mod = "_detected"; + std::string num = "_" + std::to_string(i); std::string sym = nuc.isotopicSymbol; - std::string ke_vs_th_name = sym + "_ke_vs_theta"; + std::string ke_vs_th_name = sym + num + "_ke_vs_theta"; std::string ke_vs_th_title = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)"; - std::string ke_vs_th_name_det = sym + mod + "_ke_vs_theta"; + std::string ke_vs_th_name_det = sym + num + mod + "_ke_vs_theta"; std::string ke_vs_th_title_det = ke_vs_th_name + ";#theta_{lab} (degrees);Kinetic Energy (MeV)"; - std::string ke_vs_ph_name = sym + "_ke_vs_phi"; + std::string ke_vs_ph_name = sym + num + "_ke_vs_phi"; std::string ke_vs_ph_title = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)"; - std::string ke_vs_ph_name_det = sym + mod + "_ke_vs_phi"; + std::string ke_vs_ph_name_det = sym + num + mod + "_ke_vs_phi"; std::string ke_vs_ph_title_det = ke_vs_ph_name + ";#phi_{lab} (degrees);Kinetic Energy (MeV)"; - std::string th_vs_ph_name = sym + "_theta_vs_phi"; + std::string th_vs_ph_name = sym + num + "_theta_vs_phi"; std::string th_vs_ph_title = th_vs_ph_name + ";#theta_{lab};#phi_{lab}"; - std::string th_vs_ph_name_det = sym + mod + "_theta_vs_phi"; + std::string th_vs_ph_name_det = sym + num + mod + "_theta_vs_phi"; std::string th_vs_ph_title_det = th_vs_ph_name + ";#theta_{lab};#phi_{lab}"; - std::string ex_name = sym + "_ex"; + std::string ex_name = sym + num + "_ex"; std::string ex_title = ex_name + ";E_{ex} (MeV);counts"; - std::string ex_name_det = sym + mod + "_ex"; + std::string ex_name_det = sym + num + mod + "_ex"; std::string ex_title_det = ex_name + ";E_{ex} (MeV);counts"; - std::string angdist_name = sym +"_angDist"; + std::string angdist_name = sym + num + "_angDist"; std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts"; - std::string angdist_name_det = sym + mod +"_angDist"; + std::string angdist_name_det = sym + num + mod +"_angDist"; std::string angdist_title_det = angdist_name+";cos#right(#theta_{CM}#left);counts"; - std::string hist_ke_th_name = sym + "_hist_ke_vs_theta"; + std::string hist_ke_th_name = sym + num + "_hist_ke_vs_theta"; std::string hist_ke_th_title = hist_ke_th_name + ";#theta_{lab};Kinetic Energy (MeV)"; - std::string hist_ke_th_name_det = sym + mod + "_hist_ke_vs_theta"; + std::string hist_ke_th_name_det = sym + num + mod + "_hist_ke_vs_theta"; std::string hist_ke_th_title_det = hist_ke_th_name + ";#theta_{lab};Kinetic Energy (MeV)"; MyFill(ke_vs_th_name, ke_vs_th_title, nuc.vec4.Theta()*s_rad2deg, nuc.GetKE(), 2); @@ -105,6 +114,93 @@ void RootPlotter::FillData(const Mask::Nucleus& nuc) } +//Correlation analysis for coupled step reactions, will need to be customized +//for a given reaction. By default it is off. +void RootPlotter::Correlations(const std::vector& event) +{ + static constexpr double lowerBoundEx = -2.7; + static constexpr double upperBoundEx = 2.7; + static constexpr double targetDepth = 0.5; + const Mask::Nucleus& secondary = event[6]; + const Mask::Nucleus& primary = event[4]; + const Mask::Nucleus& parent = event[3]; + const Mask::Nucleus& li = event[5]; + if(secondary.vec4.P() == 0.0 || primary.vec4.P() == 0.0) + { + return; + } + ROOT::Math::Boost boostParent(parent.vec4.BoostToCM()); + ROOT::Math::Boost boostIntermediate(li.vec4.BoostToCM()); + ROOT::Math::PxPyPzEVector a2Vec = (secondary.vec4); + ROOT::Math::PxPyPzEVector a1Vec = (primary.vec4); + + double a2KE = secondary.GetKE() + m_target.GetReverseEnergyLossFractionalDepth(secondary.Z, secondary.A, secondary.GetKE(), a2Vec.Theta(), targetDepth); + double a2P = std::sqrt(a2KE * (a2KE + 2.0 * secondary.groundStateMass)); + a2Vec.SetPxPyPzE( + a2P * std::sin(a2Vec.Theta()) * std::cos(a2Vec.Phi()), + a2P * std::sin(a2Vec.Theta()) * std::sin(a2Vec.Phi()), + a2P * std::cos(a2Vec.Theta()), + a2KE + secondary.groundStateMass + ); + + double a1KE = primary.GetKE() + m_target.GetReverseEnergyLossFractionalDepth(primary.Z, primary.A, primary.GetKE(), a1Vec.Theta(), targetDepth); + double a1P = std::sqrt(a1KE * (a1KE + 2.0 * primary.groundStateMass)); + a1Vec.SetPxPyPzE( + a1P * std::sin(a1Vec.Theta()) * std::cos(a1Vec.Phi()), + a1P * std::sin(a1Vec.Theta()) * std::sin(a1Vec.Phi()), + a1P * std::cos(a1Vec.Theta()), + a1KE + primary.groundStateMass + ); + + ROOT::Math::XYZVector a1PVec; + ROOT::Math::XYZVector a2PVec; + + ROOT::Math::PxPyPzEVector boostedIntA1Vec = boostIntermediate * a1Vec; + ROOT::Math::PxPyPzEVector boostedIntA2Vec = boostIntermediate * a2Vec; + ROOT::Math::PxPyPzEVector boostedParA2Vec = boostParent * a2Vec; + a1PVec.SetXYZ(boostedIntA1Vec.Px(), boostedIntA1Vec.Py(), boostedIntA1Vec.Pz()); + a2PVec.SetXYZ(boostedIntA2Vec.Px(), boostedIntA2Vec.Py(), boostedIntA2Vec.Pz()); + double relCosThetaCM = a1PVec.Dot(a2PVec) / (a1PVec.R() * a2PVec.R()); + + ROOT::Math::PxPyPzEVector badLi = parent.vec4 - a2Vec; + double secondaryExcitation = badLi.M() - li.groundStateMass; + + MyFill("seoncdaryEx","secondaryEx",1000,-5.0,5.0,secondaryExcitation); + MyFill("primaryEx","primaryEx",1000,-5.0,5.0,li.GetExcitationEnergy()); + MyFill("relativeDist_full","relativeDist_full",20,-1.0,1.0, relCosThetaCM); + MyFill("primaryDist_full","primaryDist_full",20,-1.0,1.0,std::cos(primary.thetaCM)); + MyFill("seondaryDist_full","seconaryDist_full",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta())); + MyFill("secondaryDist_fullOwn","seondaryDist_fullOwn",20,-1.0,1.0,std::cos(secondary.thetaCM)); + + bool isPrimary = false; + bool isSecondary = false; + if(primary.isDetected) + { + MyFill("primaryExMatched","primaryExMatched",1000,-5.0,5.0,li.GetExcitationEnergy()); + MyFill("primaryDist","primaryDist",20,-1.0,1.0,std::cos(primary.thetaCM)); + isPrimary = true; + } + if( secondary.isDetected && secondaryExcitation > lowerBoundEx && secondaryExcitation < upperBoundEx ) + { + MyFill("secondaryDist","secondaryDist",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta())); + MyFill("secondaryExMatched","secondaryExMatched",1000,-5.0,5.0,secondaryExcitation); + isSecondary = true; + } + + if(isSecondary && isPrimary) + { + return; + } + else if (isSecondary) + { + MyFill("completeDist","completeDist",20,-1.0,1.0,std::cos(boostedParA2Vec.Theta())); + } + else if (isPrimary) + { + MyFill("completeDist","completeDist",20,-1.0,1.0,std::cos(primary.thetaCM)); + } +} + void RootPlotter::MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val) { auto iter = m_map.find(name); diff --git a/src/Plotters/RootPlotter.h b/src/Plotters/RootPlotter.h index 1cc263e..19e4131 100644 --- a/src/Plotters/RootPlotter.h +++ b/src/Plotters/RootPlotter.h @@ -6,6 +6,7 @@ #include #include "Nucleus.h" +#include "Target.h" #include #include @@ -20,7 +21,8 @@ public: void Run(const std::string& inputname, const std::string& outputname); private: - void FillData(const Mask::Nucleus& nuc); + void FillData(const Mask::Nucleus& nuc, int i); + void Correlations(const std::vector& event); void MyFill(const std::string& name, const std::string& title, int bins, float min, float max, double val); void MyFill(const std::string& name, const std::string& title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy); @@ -29,6 +31,7 @@ private: std::unordered_map> m_map; static constexpr double s_rad2deg = 180.0/M_PI; + Mask::Target m_target; }; #endif \ No newline at end of file