diff --git a/include/EnergyLoss.h b/include/EnergyLoss.h index bb0bad9..969e524 100644 --- a/include/EnergyLoss.h +++ b/include/EnergyLoss.h @@ -28,7 +28,7 @@ namespace Mask { double GetEnergyLoss(int zp, int ap, double e_initial, double thickness); double GetReverseEnergyLoss(int zp, int ap, double e_final, double thickness); double GetRange(double energy); - void SetTargetComponents(std::vector& Zt, std::vector& At, std::vector& Stoich); + void SetTargetComponents(const std::vector& Zt, const std::vector& At, const std::vector& Stoich); private: double GetElectronicStoppingPower(double energy); diff --git a/include/MaskApp.h b/include/MaskApp.h index 2fb3ecd..39ec160 100644 --- a/include/MaskApp.h +++ b/include/MaskApp.h @@ -17,23 +17,19 @@ namespace Mask { bool LoadConfig(const std::string& filename); bool SaveConfig(const std::string& filename); inline int GetNumberOfSamples() const { return m_nsamples; }; - inline const std::string GetSystemName() const { return sys == nullptr ? "" : sys->GetSystemEquation(); }; + inline const std::string GetSystemName() const { return m_sys == nullptr ? "" : m_sys->GetSystemEquation(); }; inline const std::string GetOutputName() const { return m_outfile_name; }; inline const RxnType GetReactionType() const { return m_rxn_type; }; void Run(); private: - void RunOneStepRxn(); - void RunOneStepDecay(); - void RunTwoStep(); - void RunThreeStep(); - ReactionSystem* sys; + ReactionSystem* m_sys; std::string m_outfile_name; RxnType m_rxn_type; - uint32_t m_nsamples; + int m_nsamples; }; diff --git a/include/MaskFile.h b/include/MaskFile.h index afe37f9..42ffced 100644 --- a/include/MaskFile.h +++ b/include/MaskFile.h @@ -49,8 +49,8 @@ namespace Mask { FileType file_type; std::string filename; - uint32_t buffer_position; - uint32_t buffer_end; + uint64_t buffer_position; + uint64_t buffer_end; uint32_t data_size; RxnType m_rxn_type; uint32_t buffersize_bytes; diff --git a/include/RootPlotter.h b/include/RootPlotter.h index 2ea5559..17dc408 100644 --- a/include/RootPlotter.h +++ b/include/RootPlotter.h @@ -28,7 +28,7 @@ public: GenerateGraphs(); return table; }; - void FillData(const Mask::Nucleus& nuc, const std::string& modifier = ""); + void FillData(const Mask::Nucleus& nuc, double detKE = 0.0, const std::string& modifier = ""); private: THashTable* table; diff --git a/include/RxnType.h b/include/RxnType.h index 8fe3f21..bd2f8bc 100644 --- a/include/RxnType.h +++ b/include/RxnType.h @@ -36,6 +36,7 @@ namespace Mask { case RxnType::OneStepRxn: return "OneStepRxn"; case RxnType::TwoStepRxn: return "TwoStepRxn"; case RxnType::ThreeStepRxn: return "ThreeStepRxn"; + case RxnType::None : return "None"; } return "None"; @@ -51,8 +52,6 @@ namespace Mask { return static_cast(type); } - - } #endif \ No newline at end of file diff --git a/premake5.lua b/premake5.lua index 9f66766..6eba8f6 100644 --- a/premake5.lua +++ b/premake5.lua @@ -1,7 +1,7 @@ workspace "Mask" configurations { - "Debug", - "Release" + "Release", + "Debug" } project "Mask" @@ -44,8 +44,8 @@ project "RootPlot" } --User specified path to ROOT CERN libraries-- - ROOTIncludepath = "/home/gordon/cern/root-6.22.02/root-install/include" - ROOTLibpath = "/home/gordon/cern/root-6.22.02/root-install/lib" + ROOTIncludepath = "/usr/include/root" + ROOTLibpath = "/usr/lib64/root" includedirs { "include" @@ -106,4 +106,4 @@ project "DetectEff" symbols "On" filter "configurations:Release" - optimize "On" \ No newline at end of file + optimize "On" diff --git a/src/Detectors/AnasenEfficiency.cpp b/src/Detectors/AnasenEfficiency.cpp index 7998545..3f68fc7 100644 --- a/src/Detectors/AnasenEfficiency.cpp +++ b/src/Detectors/AnasenEfficiency.cpp @@ -206,48 +206,40 @@ DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus) { DetectorResult observation; //Mask::Vec3 coords; - double thetaIncident, eloss, e_dep; + double thetaIncident; for(auto& sx3 : m_Ring1) { auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi()); if(result.first != -1) { - //coords = sx3.GetHitCoordinates(result.first, result.second); observation.detectFlag = true; observation.direction = sx3.GetHitCoordinates(result.first, result.second); thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR()); if(thetaIncident > M_PI/2.0) thetaIncident = M_PI - thetaIncident; - //e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); observation.det_name = "R1"; - //return std::make_pair(true, e_dep); return observation; } } - //return std::make_pair(false, 0.0); return observation; } DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) { DetectorResult observation; - //Mask::Vec3 coords; - double thetaIncident, eloss, e_dep; + double thetaIncident; for(auto& sx3 : m_Ring2) { auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi()); if(result.first != -1) { - //coords = sx3.GetHitCoordinates(result.first, result.second); observation.detectFlag = true; observation.direction = sx3.GetHitCoordinates(result.first, result.second); thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR()); if(thetaIncident > M_PI/2.0) thetaIncident = M_PI - thetaIncident; - //e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); observation.det_name = "R2"; - //return std::make_pair(true, e_dep); return observation; } } @@ -258,22 +250,17 @@ DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) { DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) { DetectorResult observation; - //Mask::Vec3 coords; - double thetaIncident, eloss, e_dep; - + double thetaIncident; for(auto& qqq : m_forwardQQQs) { auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi()); if(result.first != -1) { - //coords = qqq.GetHitCoordinates(result.first, result.second); observation.detectFlag = true; observation.direction = qqq.GetHitCoordinates(result.first, result.second); thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR()); if(thetaIncident > M_PI/2.0) thetaIncident = M_PI - thetaIncident; - //e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); - //return std::make_pair(true, e_dep); observation.det_name = "FQQQ"; return observation; } @@ -283,22 +270,18 @@ DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) { for(auto& qqq : m_backwardQQQs) { auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi()); if(result.first != -1) { - //coords = qqq.GetHitCoordinates(result.first, result.second); observation.detectFlag = true; observation.direction = qqq.GetHitCoordinates(result.first, result.second); thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR()); if(thetaIncident > M_PI/2.0) thetaIncident = M_PI - thetaIncident; - //e_dep = det_silicon.getEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); observation.energy_deposited = det_silicon.GetEnergyLossTotal(nucleus.GetZ(), nucleus.GetA(), nucleus.GetKE(), thetaIncident); - //return std::make_pair(true, e_dep); observation.det_name = "BQQQ"; return observation; } } - //return std::make_pair(false, 0.0); return observation; } diff --git a/src/Detectors/DetectorEfficiency.cpp b/src/Detectors/DetectorEfficiency.cpp index 9680823..9d63661 100644 --- a/src/Detectors/DetectorEfficiency.cpp +++ b/src/Detectors/DetectorEfficiency.cpp @@ -25,10 +25,10 @@ int main(int argc, char** argv) { */ try { - AnasenEfficiency anasen; - anasen.CalculateEfficiency(inputname, outputname, statsname); + AnasenEfficiency anasen; + anasen.CalculateEfficiency(inputname, outputname, statsname); //std::cout<<"Running consistency check(1=success): "<& Zt, std::vector& At, std::vector& Stoich) { + void EnergyLoss::SetTargetComponents(const std::vector& Zt, const std::vector& At, const std::vector& Stoich) { comp_denom = 0; ZT = Zt; AT = At; diff --git a/src/MaskApp.cpp b/src/MaskApp.cpp index bfa090f..52e593d 100644 --- a/src/MaskApp.cpp +++ b/src/MaskApp.cpp @@ -6,14 +6,14 @@ namespace Mask { MaskApp::MaskApp() : - sys(nullptr) + m_sys(nullptr) { std::cout<<"----------GWM Kinematics Simulation----------"<>z>>a; @@ -53,7 +53,7 @@ namespace Mask { } case RxnType::OneStepRxn: { - sys = new OneStepSystem(); + m_sys = new OneStepSystem(); m_rxn_type = RxnType::OneStepRxn; for(int i=0; i<3; i++) { input>>z>>a; @@ -64,7 +64,7 @@ namespace Mask { } case RxnType::TwoStepRxn: { - sys = new TwoStepSystem(); + m_sys = new TwoStepSystem(); m_rxn_type = RxnType::TwoStepRxn; for(int i=0; i<4; i++) { input>>z>>a; @@ -75,7 +75,7 @@ namespace Mask { } case RxnType::ThreeStepRxn: { - sys = new ThreeStepSystem(); + m_sys = new ThreeStepSystem(); m_rxn_type = RxnType::TwoStepRxn; for(int i=0; i<5; i++) { input>>z>>a; @@ -87,7 +87,7 @@ namespace Mask { default: return false; } - sys->SetNuclei(zvec, avec); + m_sys->SetNuclei(zvec, avec); int nlayers; double thickness; @@ -110,7 +110,7 @@ namespace Mask { input>>z>>a>>s; zvec.push_back(z); avec.push_back(a); svec.push_back(s); } - sys->AddTargetLayer(zvec, avec, svec, thickness); + m_sys->AddTargetLayer(zvec, avec, svec, thickness); input>>junk; } std::cout<<"Reaction equation: "<>junk>>m_nsamples; input>>junk>>par1>>junk>>par2; - sys->SetBeamDistro(par1, par2); + m_sys->SetBeamDistro(par1, par2); input>>junk>>par1; switch(m_rxn_type) { case RxnType::OneStepRxn : { - dynamic_cast(sys)->SetReactionThetaType(par1); + dynamic_cast(m_sys)->SetReactionThetaType(par1); break; } case RxnType::TwoStepRxn : { - dynamic_cast(sys)->SetReactionThetaType(par1); + dynamic_cast(m_sys)->SetReactionThetaType(par1); break; } case RxnType::ThreeStepRxn : { - dynamic_cast(sys)->SetReactionThetaType(par1); + dynamic_cast(m_sys)->SetReactionThetaType(par1); break; } } input>>junk>>par1>>junk>>par2; - sys->SetTheta1Range(par1, par2); + m_sys->SetTheta1Range(par1, par2); input>>junk>>par1>>junk>>par2; - sys->SetPhi1Range(par1, par2); + m_sys->SetPhi1Range(par1, par2); input>>junk>>par1>>junk>>par2; - sys->SetExcitationDistro(par1, par2); + m_sys->SetExcitationDistro(par1, par2); input>>junk>>dfile1; input>>junk>>dfile2; switch(m_rxn_type) { case RxnType::OneStepRxn : { - DecaySystem* this_sys = dynamic_cast(sys); + DecaySystem* this_sys = dynamic_cast(m_sys); this_sys->SetDecay1Distribution(dfile1); std::cout<<"Decay1 angular momentum: "<GetDecay1AngularMomentum()<GetDecay1BranchingRatio()<(sys); + TwoStepSystem* this_sys = dynamic_cast(m_sys); this_sys->SetDecay1Distribution(dfile1); std::cout<<"Decay1 angular momentum: "<GetDecay1AngularMomentum()<GetDecay1BranchingRatio()<(sys); + ThreeStepSystem* this_sys = dynamic_cast(m_sys); this_sys->SetDecay1Distribution(dfile1); this_sys->SetDecay2Distribution(dfile2); std::cout<<"Decay1 angular momentum: "<GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<GetDecay2AngularMomentum()<RunSystem(); - output.WriteData(sys->GetNuclei()); + m_sys->RunSystem(); + output.WriteData(m_sys->GetNuclei()); } output.Close(); diff --git a/src/MaskFile.cpp b/src/MaskFile.cpp index 8c95eae..23de319 100644 --- a/src/MaskFile.cpp +++ b/src/MaskFile.cpp @@ -332,7 +332,7 @@ namespace Mask { } } - unsigned int local_end = buffer_position + data_size; + uint64_t local_end = buffer_position + data_size; if(local_end > buffer_end) { std::cerr<<"Attempting to read past end of file at MaskFile::ReadData! Returning empty"<