From 87aa3dba8a3f195d32cf7b0cf745412e9f6d1f3d Mon Sep 17 00:00:00 2001 From: Gordon McCann Date: Mon, 23 Aug 2021 10:00:36 -0400 Subject: [PATCH] Major overhaul. Fixed bugs in decay angular distribution, now classed. Added ANASEN geometry. Minor bugfixes throughout. Added center of mass switch for main reaction (better randomization where ejectile is not fixed in theta). --- etc/12C3Hed_pdecay_L2_fit.txt | 4 + etc/12C3Hed_pdecay_L2_theory.txt | 4 + etc/isotropic_dist.txt | 2 + include/AnasenEfficiency.h | 60 ++++ include/AngularDistribution.h | 37 ++ include/KinematicsExceptions.h | 22 +- include/LegendrePoly.h | 6 + include/MaskFile.h | 54 +++ include/QQQDetector.h | 48 +++ include/Reaction.h | 13 +- include/ReactionSystem.h | 21 +- include/SabreEfficiency.h | 21 +- include/StripDetector.h | 60 ++++ input.txt | 21 +- src/AnasenEfficiency.cpp | 562 +++++++++++++++++++++++++++++++ src/AngularDistribution.cpp | 104 ++++++ src/EnergyLoss.cpp | 11 +- src/Kinematics.cpp | 34 +- src/LegendrePoly.cpp | 23 +- src/MaskFile.cpp | 179 ++++++++++ src/Plotter.cpp | 4 +- src/QQQDetector.cpp | 155 +++++++++ src/Reaction.cpp | 86 ++++- src/ReactionSystem.cpp | 37 +- src/SabreEfficiency.cpp | 487 +++++++++++--------------- src/StripDetector.cpp | 141 ++++++++ src/ThreeStepSystem.cpp | 21 +- src/TwoStepSystem.cpp | 10 +- src/main.cpp | 57 +--- 29 files changed, 1890 insertions(+), 394 deletions(-) create mode 100644 etc/12C3Hed_pdecay_L2_fit.txt create mode 100644 etc/12C3Hed_pdecay_L2_theory.txt create mode 100644 etc/isotropic_dist.txt create mode 100644 include/AnasenEfficiency.h create mode 100644 include/AngularDistribution.h create mode 100644 include/MaskFile.h create mode 100644 include/QQQDetector.h create mode 100644 include/StripDetector.h create mode 100644 src/AnasenEfficiency.cpp create mode 100644 src/AngularDistribution.cpp create mode 100644 src/MaskFile.cpp create mode 100644 src/QQQDetector.cpp create mode 100644 src/StripDetector.cpp diff --git a/etc/12C3Hed_pdecay_L2_fit.txt b/etc/12C3Hed_pdecay_L2_fit.txt new file mode 100644 index 0000000..972b20c --- /dev/null +++ b/etc/12C3Hed_pdecay_L2_fit.txt @@ -0,0 +1,4 @@ +AngularMomentum: 2 +A0: 0.2121 +A2: 0.2910 +A4: 0.03636 diff --git a/etc/12C3Hed_pdecay_L2_theory.txt b/etc/12C3Hed_pdecay_L2_theory.txt new file mode 100644 index 0000000..c3543a0 --- /dev/null +++ b/etc/12C3Hed_pdecay_L2_theory.txt @@ -0,0 +1,4 @@ +AngularMomentum: 2 +A0: 0.5 +A2: 0.544 +A4: 0.362 diff --git a/etc/isotropic_dist.txt b/etc/isotropic_dist.txt new file mode 100644 index 0000000..2b1fe94 --- /dev/null +++ b/etc/isotropic_dist.txt @@ -0,0 +1,2 @@ +AngularMomentum: 0 +A0: 0.5 \ No newline at end of file diff --git a/include/AnasenEfficiency.h b/include/AnasenEfficiency.h new file mode 100644 index 0000000..054a99c --- /dev/null +++ b/include/AnasenEfficiency.h @@ -0,0 +1,60 @@ +#ifndef ANASEN_EFFICIENCY_H +#define ANASEN_EFFICIENCY_H + +#include +#include + +#include "StripDetector.h" +#include "QQQDetector.h" + +class AnasenEfficiency { +public: + AnasenEfficiency(); + ~AnasenEfficiency(); + void CalculateEfficiency(const std::string& filename); + inline void SetReactionType(int type) { m_rxn_type = type; }; + void DrawDetectorSystem(const std::string& filename); + double RunConsistencyCheck(); + +private: + void MyFill(THashTable* table, const std::string& name, const std::string& title, int bins, float min, float max, double val); + void MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy); + void RunDecay(const std::string& filename); + void RunTwoStep(const std::string& filename); + void RunThreeStep(const std::string& filename); + + bool IsRing1(double theta, double phi); + bool IsRing2(double theta, double phi); + bool IsQQQ(double theta, double phi); + + inline bool IsDoubleEqual(double x, double y) { return std::fabs(x-y) < epsilon ? true : false; }; + + int m_rxn_type; + std::vector m_Ring1, m_Ring2; + std::vector m_forwardQQQs; + std::vector m_backwardQQQs; + + /**** ANASEN geometry constants *****/ + const int n_sx3_per_ring = 12; + const int n_qqq = 4; + const double sx3_length = 0.075; + const double sx3_width = 0.04; + const double barrel_gap = 0.013 + 0.049; //0.049 is base gap due to frames + const double ring1_z = sx3_length/2.0 + barrel_gap/2.0; + //const double ring2_z = -0.124 + sx3_length/2.0 + 0.0245 - barrel_gap/2.0; + const double qqq_nom_z = 0.025 + sx3_length + 0.0245 + barrel_gap/2.0; + const double qqq_rinner = 0.0501; + const double qqq_router = 0.0990; + const double qqq_deltaphi = 1.52119; + const double qqq_z[4] = {qqq_nom_z, qqq_nom_z - 0.00828, qqq_nom_z, qqq_nom_z}; + const double qqq_phi[4] = {5.49779, 0.785398, 2.35619, 3.92699}; + const double ring_rho[12] = {0.0890601, 0.0889871, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247, 0.0890354, 0.0890354}; + const double ring_phi[12] = {0.785795, 0.262014, 6.02132, 5.49779, 4.97426, 4.45052, 3.92699, 3.40346, 2.87972, 2.35619, 1.83266, 1.30893}; + /*************************/ + + static constexpr double threshold = 0.2; //MeV + static constexpr double epsilon = 1.0e-6; + static constexpr double deg2rad = M_PI/180.0; +}; + +#endif \ No newline at end of file diff --git a/include/AngularDistribution.h b/include/AngularDistribution.h new file mode 100644 index 0000000..ef9f3a8 --- /dev/null +++ b/include/AngularDistribution.h @@ -0,0 +1,37 @@ +#ifndef ANGULARDISTRIBUTION_H +#define ANGULARDISTRIBUTION_H + +#include +#include +#include + +class AngularDistribution { +public: + AngularDistribution(); + AngularDistribution(const std::string& file); + ~AngularDistribution(); + void ReadDistributionFile(const std::string& file); + void AttachRandomNumberGenerator(TRandom3* random) { generator = random; }; + double GetRandomCosTheta(); + int GetL() { return L; }; + double GetBranchingRatio() { return branchingRatio; }; + +private: + bool IsIsotropic(); + bool IsGeneratorSet() { + if(generator) { + return true; + } else { + return false; + } + } + + TRandom3* generator; //NOT OWNED BY ANGULAR DISTRIBUTION + + double branchingRatio; + int L; + std::vector constants; + bool isoFlag; +}; + +#endif \ No newline at end of file diff --git a/include/KinematicsExceptions.h b/include/KinematicsExceptions.h index e59972f..76ed6f1 100644 --- a/include/KinematicsExceptions.h +++ b/include/KinematicsExceptions.h @@ -2,6 +2,7 @@ #define KINEMATICSEXCEPTIONS_H #include +#include #include /* ELossException @@ -12,8 +13,16 @@ for locations where this could be thrown */ struct ELossException : public std::exception { + ELossException(const std::string& error) { + m_error = error; + } + + std::string m_error; + const char* what() const noexcept { - return "Failure to calculate particle energy loss. See KinematicsExceptions.h for documentation."; + std::string err_str = "Failure to calculate particle energy loss. Reason: "; + err_str += m_error + " See KinematicsExceptions.h for documentation."; + return err_str.c_str(); }; }; @@ -61,6 +70,17 @@ struct QValueException : public std::exception { }; }; +/* + EnergyThresholdException + This is an exception thrown when the Reaction attempts to calculate a reaction which does not have enough incoming (read: beam kinetic energy) energy + to occur even for the 0 degree case. The reaction is kinematically forbidden. +*/ +struct EnergyThresholdException : public std::exception { + const char* what() const noexcept { + return "Reaction does not have enough energy to proceed. See KinematicsExceptions.h for documentation."; + }; +}; + #endif \ No newline at end of file diff --git a/include/LegendrePoly.h b/include/LegendrePoly.h index d05808b..270d42e 100644 --- a/include/LegendrePoly.h +++ b/include/LegendrePoly.h @@ -2,5 +2,11 @@ #define LEGENDREPOLY_H double P_l(int l, double x); +double P_l_ROOT(double* x, double* pars); +double Normed_P_l_sq(int l, double x); + +double P_0(double x); +double P_1(double x); +double P_2(double x); #endif \ No newline at end of file diff --git a/include/MaskFile.h b/include/MaskFile.h new file mode 100644 index 0000000..b584f89 --- /dev/null +++ b/include/MaskFile.h @@ -0,0 +1,54 @@ +#ifndef MASKFILE_H +#define MASKFILE_H + +#include +#include +#include + +#include "Nucleus.h" + +namespace Mask { + +class MaskFile { +public: + enum class FileType { + read, + write, + append, + none + }; + + MaskFile(); + MaskFile(const std::string& name, MaskFile::FileType type); + bool Open(const std::string& name, MaskFile::FileType type); + inline bool IsOpen() { return file.is_open(); } + void Close(); + + void WriteHeader(std::vector& data); + void WriteData(std::vector& data); + void ReadHeader(); + std::vector ReadData(); + + + +private: + void FlushBuffer(); + void FillBuffer(); + + FileType file_type; + std::string filename; + unsigned int buffer_position; + int count; + + std::fstream file; + + std::vector> data_buffer; + + static constexpr int buffersize = 1000; + static constexpr int width = 0; + static constexpr int precision = 3; +}; + +}; + +#endif \ No newline at end of file diff --git a/include/QQQDetector.h b/include/QQQDetector.h new file mode 100644 index 0000000..3ffee56 --- /dev/null +++ b/include/QQQDetector.h @@ -0,0 +1,48 @@ +/* + QQQDetector.h + Class implementing geometry for QQQDetector where the detector is perpendicular to the beam axis. + Detector is first generated centered on the x-axis (phi=0) + + Coordinate convention : +z is downstream, -z is upstream. +y is vertically up in the lab. +*/ +#ifndef QQQDETECTOR_H +#define QQQDETECTOR_H + +#include +#include + +#include "Vec3.h" +#include "Rotation.h" + +class QQQDetector { +public: + QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x=0, double y=0); + ~QQQDetector(); + inline Mask::Vec3 GetRingCoordinates(int ringch, int corner) { return m_ringCoords[ringch][corner]; }; + inline Mask::Vec3 GetWedgeCoordinates(int wedgech, int corner) { return m_wedgeCoords[wedgech][corner]; }; + Mask::Vec3 GetTrajectoryCoordinates(double theta, double phi); + std::pair GetTrajectoryRingWedge(double theta, double phi); + Mask::Vec3 GetHitCoordinates(int ringch, int wedgech); + + inline int GetNumberOfRings() { return nrings; }; + inline int GetNumberOfWedges() { return nwedges; }; + +private: + + inline bool CheckChannel(int ch) { return (ch >=0 && ch < nrings); }; + inline bool CheckCorner(int corner) { return (corner >=0 && corner < 4); }; + + void CalculateCorners(); + Mask::Vec3 TransformCoordinates(Mask::Vec3& vector) { return m_ZRot*vector + m_translation; }; + + double m_Rinner, m_Router, m_deltaR, m_deltaPhi, m_deltaPhi_per_wedge, m_phiCentral; + std::vector> m_ringCoords, m_wedgeCoords; + Mask::Vec3 m_translation; + Mask::ZRotation m_ZRot; + + static constexpr int nrings = 16; + static constexpr int nwedges = 16; + static constexpr double deg2rad = M_PI/180.0; +}; + +#endif \ No newline at end of file diff --git a/include/Reaction.h b/include/Reaction.h index d8b1b40..1c3b940 100644 --- a/include/Reaction.h +++ b/include/Reaction.h @@ -23,6 +23,7 @@ public: void SetNuclei(int zt, int at, int zp, int ap, int ze, int ae); void SetNuclei(const Nucleus* nucs); void SetBeamKE(double bke); + void SetEjectileThetaType(int type); /*Setters and getters*/ inline void SetLayeredTarget(LayeredTarget* targ) { target = targ; }; @@ -46,10 +47,16 @@ public: inline const Nucleus& GetTarget() const { return reactants[0]; }; inline const Nucleus& GetEjectile() const { return reactants[2]; }; inline const Nucleus& GetResidual() const { return reactants[3]; }; + inline void ResetTarget() { reactants[0].SetVectorCartesian(0,0,0,0); }; + inline void ResetProjectile() { reactants[1].SetVectorCartesian(0,0,0,0); }; + inline void ResetEjectile() { reactants[2].SetVectorCartesian(0,0,0,0); }; + inline void ResetResidual() { reactants[3].SetVectorCartesian(0,0,0,0); }; inline int GetRxnLayer() { return rxnLayer; }; private: void CalculateReaction(); //target + project -> eject + resid + void CalculateReactionThetaLab(); + void CalculateReactionThetaCM(); void CalculateDecay(); //target -> light_decay (eject) + heavy_decay(resid) Nucleus reactants[4]; //0=target, 1=projectile, 2=ejectile, 3=residual @@ -57,10 +64,14 @@ private: double m_bke, m_theta, m_phi, m_ex; - int rxnLayer; + int rxnLayer; + int m_eject_theta_type; bool decayFlag, nuc_initFlag, resid_elossFlag; + static constexpr int lab = 0; + static constexpr int center_of_mass = 1; + }; }; diff --git a/include/ReactionSystem.h b/include/ReactionSystem.h index dde2105..10d30dc 100644 --- a/include/ReactionSystem.h +++ b/include/ReactionSystem.h @@ -10,6 +10,7 @@ #define REACTIONSYSTEM_H #include "Reaction.h" +#include "AngularDistribution.h" #include #include @@ -25,15 +26,14 @@ public: 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; }; + void SetRandomGenerator(TRandom3* gen); inline void SetBeamDistro(double mean, double sigma) { m_beamDist = std::make_pair(mean, sigma); }; + inline void SetReactionThetaType(int type) { step1.SetEjectileThetaType(type); }; inline void SetTheta1Range(double min, double max) { m_theta1Range = std::make_pair(min*deg2rad, max*deg2rad); }; + inline void SetPhi1Range(double min, double max) { m_phi1Range = std::make_pair(min*deg2rad, max*deg2rad); }; inline void SetExcitationDistro(double mean, double sigma) { m_exDist = std::make_pair(mean, sigma); }; - inline void SetDecay1AngularMomentum(double l) { L1 = l; }; - inline void SetDecay2AngularMomentum(double l) { L2 = l; }; - - /*Sampling over angular distribution*/ - double GetDecayTheta(int L); + inline void SetDecay1Distribution(const std::string& file) { decay1dist.ReadDistributionFile(file); }; + inline void SetDecay2Distribution (const std::string& file) { decay2dist.ReadDistributionFile(file); }; virtual void RunSystem(); @@ -42,6 +42,10 @@ public: inline const Nucleus& GetEjectile() const { return step1.GetEjectile(); }; inline const Nucleus& GetResidual() const { return step1.GetResidual(); }; inline const char* GetSystemEquation() const { return m_sys_equation.c_str(); }; + inline int GetDecay1AngularMomentum() { return decay1dist.GetL(); }; + inline int GetDecay2AngularMomentum(){ return decay2dist.GetL(); }; + inline double GetDecay1BranchingRatio() { return decay1dist.GetBranchingRatio(); }; + inline double GetDecay2BranchingRatio(){ return decay2dist.GetBranchingRatio(); }; protected: virtual void LinkTarget(); @@ -51,12 +55,13 @@ protected: LayeredTarget target; //Sampling information - std::pair m_beamDist, m_theta1Range, m_exDist; + std::pair m_beamDist, m_theta1Range, m_phi1Range, m_exDist; TRandom3* generator; //not owned by ReactionSystem + AngularDistribution decay1dist, decay2dist; + bool target_set_flag, gen_set_flag; int rxnLayer; - int L1, L2; std::string m_sys_equation; static constexpr double deg2rad = M_PI/180.0; }; diff --git a/include/SabreEfficiency.h b/include/SabreEfficiency.h index 67ff901..a57dd6c 100644 --- a/include/SabreEfficiency.h +++ b/include/SabreEfficiency.h @@ -4,6 +4,7 @@ #include "SabreDetector.h" #include "Target.h" #include "DeadChannelMap.h" +#include "Kinematics.h" #include class SabreEfficiency { @@ -12,20 +13,23 @@ public: ~SabreEfficiency(); inline void SetReactionType(int t) { m_rxn_type = t; }; void SetDeadChannelMap(std::string& filename) { dmap.LoadMapfile(filename); }; - void CalculateEfficiency(const char* file); + void CalculateEfficiency(const std::string& file); + void DrawDetectorSystem(const std::string& filename); + double RunConsistencyCheck(); private: - void MyFill(THashTable* table, const char* name, const char* title, int bins, float min, float max, double val); - void MyFill(THashTable* table, const char* name, const char* title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy); - void Run2Step(const char*); - void Run3Step(const char*); - void RunDecay(const char*); + void MyFill(THashTable* table, const std::string& name, const std::string& title, int bins, float min, float max, double val); + void MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy); + std::pair IsSabre(Mask::NucData* nucleus); + void Run2Step(const std::string& filename); + void Run3Step(const std::string& filename); + void RunDecay(const std::string& filename); int m_rxn_type; std::vector detectors; - std::vector ringxs, ringys, ringzs; - std::vector wedgexs, wedgeys, wedgezs; + Target deadlayer; + Target sabre_eloss; DeadChannelMap dmap; @@ -42,6 +46,7 @@ private: const double PHI4 = 90.0; const double DEG2RAD = M_PI/180.0; static constexpr double DEADLAYER_THIN = 50 * 1e-7 * 2.3296 * 1e6; // ug/cm^2 (50 nm thick * density) + static constexpr double SABRE_THICKNESS = 500 * 1e-4 * 2.3926 * 1e6; // ug/cm^2 (500 um thick * density) const double ENERGY_THRESHOLD = 0.2; //in MeV diff --git a/include/StripDetector.h b/include/StripDetector.h new file mode 100644 index 0000000..2aaef13 --- /dev/null +++ b/include/StripDetector.h @@ -0,0 +1,60 @@ +#ifndef __STRIP_DETECTOR_H +#define __STRIP_DETECTOR_H + +// +z is along beam axis +// +y is vertically "upward" in the lab frame + +//angles must be in radians, but distances can be whatever +//PROVIDED all input distances are the same + +//Front strips from lowest y to highest y + +//Back strips from lowest z to highest z + +#include +#include +#include + +#include "Vec3.h" +#include "Rotation.h" + +class StripDetector { +public: + + StripDetector(int ns, double len, double wid, double cphi, double cz, double crho); + ~StripDetector(); + inline Mask::Vec3 GetFrontStripCoordinates(int stripch, int corner) { return front_strip_coords[stripch][corner]; }; + inline Mask::Vec3 GetBackStripCoordinates(int stripch, int corner) { return back_strip_coords[stripch][corner]; }; + inline Mask::Vec3 GetRotatedFrontStripCoordinates(int stripch, int corner) { return rotated_front_strip_coords[stripch][corner]; }; + inline Mask::Vec3 GetRotatedBackStripCoordinates(int stripch, int corner) { return rotated_back_strip_coords[stripch][corner]; }; + inline void SetRandomNumberGenerator(TRandom3* random) { m_random = random; }; + Mask::Vec3 GetHitCoordinates(int front_stripch, double front_strip_ratio); + std::pair GetChannelRatio(double theta, double phi); + +private: + void CalculateCorners(); + + int num_strips; + static constexpr int num_corners = 4; + + double length; //common to all strips, hence total + double total_width; + double front_strip_width; //assuming equal widths + double back_strip_length; //assuming equal widths + + double center_phi; //assuming det centered above x-axis (corresponds to zero phi) + double center_z; + double center_rho; //perpendicular radius from axis + + std::vector> front_strip_coords, back_strip_coords; + std::vector> rotated_front_strip_coords, rotated_back_strip_coords; + + Mask::ZRotation zRot; + + TRandom3* m_random; //Not owned by StripDetector! + + inline bool ValidChannel(int f) { return ((f >= 0 && f < num_strips) ? true : false); }; + inline bool ValidRatio(double r) { return ((r >= -1 && r <= 1) ? true : false); }; +}; + +#endif diff --git a/input.txt b/input.txt index 3b38af1..768161a 100644 --- a/input.txt +++ b/input.txt @@ -1,28 +1,31 @@ ----------Data Information---------- -OutputFile: /data1/gwm17/10B3He/Feb2021/12C3Hed_13N3546keV_ps_fixedThetaBins_test.root +OutputFile: /data1/gwm17/TRIUMF_7Bed/simulation/7Bedp_600keV_beam_centered_target_targetgap_BackQQQ_rndmCM_test.root SaveTree: yes SavePlots: yes ----------Reaction Information---------- ReactionType: 2 Z A (order is target, projectile, ejectile, break1, break3) -6 12 -2 3 1 2 +4 7 1 1 +2 4 ----------Target Information---------- Name: test_targ Layers: 1 ~Layer1 -Thickness(ug/cm^2): 41 +Thickness(ug/cm^2): 50 Z A Stoich +1 2 2 6 12 1 0 ~ ----------Sampling Information---------- NumberOfSamples: 100000 -BeamMeanEnergy(MeV): 24 BeamEnergySigma(MeV): 0.001 -EjectileThetaMin(deg): 3.0 EjectileThetaMax(deg): 3.0 -ResidualExMean(MeV): 3.546 ResidualExSigma(MeV): 0.02 -Decay1_AngularMomentum: 0 -Decay2_AngularMomentum: 0 +BeamMeanEnergy(MeV): 0.6 BeamEnergySigma(MeV): 0.001 +EjectileThetaType(0=Lab,1=CM): 1 +EjectileThetaMin(deg): 0.0 EjectileThetaMax(deg): 180.0 +EjectilePhiMin(deg): 0.0 EjectilePhiMax(deg): 360.0 +ResidualExMean(MeV): 0.0 ResidualExSigma(MeV): 0.0 +Decay1_DistributionFile: ./etc/isotropic_dist.txt +Decay2_DistributionFile: ./etc/isotropic_dist.txt -------------------------------------- diff --git a/src/AnasenEfficiency.cpp b/src/AnasenEfficiency.cpp new file mode 100644 index 0000000..7367cc7 --- /dev/null +++ b/src/AnasenEfficiency.cpp @@ -0,0 +1,562 @@ +#include "AnasenEfficiency.h" +#include "Kinematics.h" +#include +#include +#include +#include +#include +#include + +AnasenEfficiency::AnasenEfficiency() { + for(int i=0; iFindObject(name.c_str()); + if(h) { + h->Fill(val); + } else { + h = new TH1F(name.c_str(), title.c_str(), bins, min, max); + h->Fill(val); + table->Add(h); + } +} + +void AnasenEfficiency::MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy) { + TH2F* h = (TH2F*) table->FindObject(name.c_str()); + if(h) { + h->Fill(valx, valy); + } else { + h = new TH2F(name.c_str(), title.c_str(), binsx, minx, maxx, binsy, miny, maxy); + h->Fill(valx, valy); + table->Add(h); + } +} + +void AnasenEfficiency::DrawDetectorSystem(const std::string& filename) { + TFile* file = TFile::Open(filename.c_str(), "RECREATE"); + + std::vector x, y, z; + std::vector cx, cy, cz; + Mask::Vec3 coords; + for(int i=0; iSetName("CornerGraph"); + graph->SetMarkerStyle(2); + graph->SetLineColor(1); + + TGraph2D* graph2 = new TGraph2D(cx.size(), &(cx[0]), &(cy[0]), &(cz[0])); + graph2->SetName("CenterGraph"); + graph2->SetMarkerStyle(2); + graph2->SetMarkerColor(4); + + TCanvas* canvas = new TCanvas(); + canvas->SetName("ANASEN Detector"); + graph->Draw("A|P"); + graph2->Draw("same|P"); + canvas->Write(); + graph->Write(); + graph2->Write(); + file->Close(); +} + +double AnasenEfficiency::RunConsistencyCheck() { + std::vector r1_points; + std::vector r2_points; + std::vector fqqq_points; + std::vector bqqq_points; + for(int i=0; iGet("DataTree"); + + THashTable* table = new THashTable(); + Mask::NucData* eject = new Mask::NucData(); + Mask::NucData* resid = new Mask::NucData(); + + tree->SetBranchAddress("ejectile", &eject); + tree->SetBranchAddress("residual", &resid); + + double nevents = tree->GetEntries(); + + //Progress tracking + int percent5 = nevents*0.05; + int count = 0; + int npercent = 0; + + int detected_eject=0; + int detected_resid=0; + for(int i=0; iGetEntries(); i++) { + if(++count == percent5) {//Show progress every 5% + npercent++; + count = 0; + std::cout<<"\rPercent completed: "<GetEntry(i); + + if(eject->KE >= threshold && (IsRing1(eject->theta, eject->phi) || IsRing2(eject->theta, eject->phi) || IsQQQ(eject->theta, eject->phi))) { + detected_eject++; + MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/deg2rad,300,0.0,30.0,eject->KE); + } + + if(resid->KE > threshold && (IsRing1(resid->theta, resid->phi) || IsRing2(resid->theta, resid->phi) || IsQQQ(resid->theta, resid->phi))) { + detected_resid++; + MyFill(table, "detected_resid_theta_ke","detected residual theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,resid->theta/deg2rad,300,0.0,30.0,resid->KE); + } + + } + + double ejecteff = ((double) detected_eject)/nevents; + double resideff = ((double) detected_resid)/nevents; + TParameter eject_eff("Light Breakup Efficiency", ejecteff); + TParameter resid_eff("Heavy Breakup Efficiency", resideff); + + input->cd(); + eject_eff.Write(); + resid_eff.Write(); + table->Write(); + input->Close(); +} + +void AnasenEfficiency::RunTwoStep(const std::string& filename) { + + TFile* input = TFile::Open(filename.c_str(), "UPDATE"); + TTree* tree = (TTree*) input->Get("DataTree"); + + Mask::NucData* eject = new Mask::NucData(); + Mask::NucData* break1 = new Mask::NucData(); + Mask::NucData* break2 = new Mask::NucData(); + + THashTable* table = new THashTable(); + + tree->SetBranchAddress("ejectile", &eject); + tree->SetBranchAddress("breakup1", &break1); + tree->SetBranchAddress("breakup2", &break2); + + double nevents = tree->GetEntries(); + + //Progress tracking + int percent5 = nevents*0.05; + int count = 0; + int npercent = 0; + + double costheta_cm =0; + + bool break1_flag, break2_flag, eject_flag; + int b1_count=0, b2_count=0, eject_count=0; + int b1b2_count=0, b1e_count=0, b2e_count=0, b1b2e_count=0; + for(int i=0; iGetEntries(); i++) { + if(++count == percent5) {//Show progress every 5% + npercent++; + count = 0; + std::cout<<"\rPercent completed: "<GetEntry(i); + + if(eject->KE >= threshold && (IsRing1(eject->theta, eject->phi) || IsRing2(eject->theta, eject->phi) || IsQQQ(eject->theta, eject->phi))) { + eject_count++; + eject_flag = true; + MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/deg2rad,300,0.0,30.0,eject->KE); + } + + if(break1->KE >= threshold && (IsRing1(break1->theta, break1->phi) || IsRing2(break1->theta, break1->phi) || IsQQQ(break1->theta, break1->phi))) { + b1_count++; + break1_flag = true; + costheta_cm = std::cos(break1->theta_cm); + MyFill(table,"detected_break1_cm_theta","detected breakup1 CM theta;cos(#theta_{CM})",20,-1.0,1.0,costheta_cm); + MyFill(table, "detected_break1_theta_ke","detected breakup1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/deg2rad,300,0.0,30.0,break1->KE); + } + + if(break2->KE >= threshold && (IsRing1(break2->theta, break2->phi) || IsRing2(break2->theta, break2->phi) || IsQQQ(break2->theta, break2->phi))) { + b2_count++; + break2_flag = true; + MyFill(table, "detected_break2_theta_ke","detected breakup2 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break2->theta/deg2rad,300,0.0,30.0,break2->KE); + } + + if(break1_flag && break2_flag && eject_flag) { + b1b2e_count++; + b1b2_count++; + b1e_count++; + b2e_count++; + MyFill(table,"b1b2_theta_theta","b1b2_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break1->theta/deg2rad,180,0.0,180.0,break2->theta/deg2rad); + MyFill(table,"b1b2_ke_ke","b1b2_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break1->KE,300,0.0,30.0,break2->KE); + MyFill(table,"b1b2_phi_phi","b1b2_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break1->phi/deg2rad,360,0.0,360.0,break2->phi/deg2rad); + } else if(break1_flag && eject_flag) { + b1e_count++; + } else if(break2_flag && eject_flag) { + b2e_count++; + } else if(break1_flag && break2_flag) { + b1b2_count++; + MyFill(table,"b1b2_theta_theta","b1b2_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break1->theta/deg2rad,180,0.0,180.0,break2->theta/deg2rad); + MyFill(table,"b1b2_ke_ke","b1b2_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break1->KE,300,0.0,30.0,break2->KE); + MyFill(table,"b1b2_phi_phi","b1b2_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break1->phi/deg2rad,360,0.0,360.0,break2->phi/deg2rad); + } + } + + double eeff = ((double) eject_count)/nevents; + double b1eff = ((double) b1_count)/nevents; + double b2eff = ((double) b2_count)/nevents; + double b1b2eff = b1b2_count/nevents; + double b1eeff = b1e_count/nevents; + double b2eeff = b2e_count/nevents; + double b1b2eeff = b1b2e_count/nevents; + TParameter eject_eff("Ejectile Efficiency", eeff); + TParameter break1_eff("Breakup 1 Efficiency", b1eff); + TParameter break2_eff("Breakup 2 Efficiency", b2eff); + TParameter break1break2_eff("Breakup1 Coincident with Breakup2", b1b2eff); + TParameter break1eject_eff("Breakup1 Coincident with Ejectile", b1eeff); + TParameter break2eject_eff("Breakup2 Coincident with Ejectile", b2eeff); + TParameter break1break2eject_eff("All Three Particles Coincident", b1b2eeff); + + + input->cd(); + table->Write(); + eject_eff.Write(); + break1_eff.Write(); + break2_eff.Write(); + break1break2_eff.Write(); + break1eject_eff.Write(); + break2eject_eff.Write(); + break1break2eject_eff.Write(); + + input->Close(); +} + +void AnasenEfficiency::RunThreeStep(const std::string& filename) { + TFile* input = TFile::Open(filename.c_str(), "UPDATE"); + TTree* tree = (TTree*) input->Get("DataTree"); + + THashTable* table = new THashTable(); + + Mask::NucData* break1 = new Mask::NucData(); + Mask::NucData* break3 = new Mask::NucData(); + Mask::NucData* break4 = new Mask::NucData(); + + + tree->SetBranchAddress("breakup1", &break1); + tree->SetBranchAddress("breakup3", &break3); + tree->SetBranchAddress("breakup4", &break4); + + double nevents = tree->GetEntries(); + + //Progress tracking + int percent5 = nevents*0.05; + int count = 0; + int npercent = 0; + + bool break1_flag, break3_flag, break4_flag; + int b1_count=0, b3_count=0, b4_count=0; + int b1b3_count=0, b1b4_count=0, b3b4_count=0, b1b3b4_count=0; + + for(int i=0; iGetEntries(); i++) { + if(++count == percent5) {//Show progress every 5% + npercent++; + count = 0; + std::cout<<"\rPercent completed: "<GetEntry(i); + + if(break1->KE >= threshold && (IsRing1(break1->theta, break1->phi) || IsRing2(break1->theta, break1->phi) || IsQQQ(break1->theta, break1->phi))) { + b1_count++; + break1_flag = true; + MyFill(table, "detected_break1_theta_ke","detected breakup1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/deg2rad,300,0.0,30.0,break1->KE); + } + + + if(break3->KE >= threshold && (IsRing1(break3->theta, break3->phi) || IsRing2(break3->theta, break3->phi) || IsQQQ(break3->theta, break3->phi))) { + b3_count++; + break3_flag = true; + MyFill(table, "detected_break3_theta_ke","detected breakup3 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break3->theta/deg2rad,300,0.0,30.0,break3->KE); + } + + if(break4->KE >= threshold && (IsRing1(break4->theta, break4->phi) || IsRing2(break4->theta, break4->phi) || IsQQQ(break4->theta, break4->phi))) { + b4_count++; + break4_flag = true; + MyFill(table, "detected_break4_theta_ke","detected breakup4 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break4->theta/deg2rad,300,0.0,30.0,break4->KE); + } + + if(break1_flag && break3_flag && break4_flag) { + b1b3b4_count++; + b1b3_count++; + b1b4_count++; + b3b4_count++; + MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/deg2rad,180,0.0,180.0,break4->theta/deg2rad); + MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break3->KE,300,0.0,30.0,break4->KE); + MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/deg2rad,360,0.0,360.0,break4->phi/deg2rad); + } else if(break1_flag && break3_flag) { + b1b3_count++; + } else if(break1_flag && break4_flag) { + b1b4_count++; + } else if(break3_flag && break4_flag) { + b3b4_count++; + MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/deg2rad,180,0.0,180.0,break4->theta/deg2rad); + MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",300,0.0,30.0,break3->KE,300,0.0,30.0,break4->KE); + MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/deg2rad,360,0.0,360.0,break4->phi/deg2rad); + } + } + + double b1eff = ((double) b1_count)/nevents; + double b3eff = ((double) b3_count)/nevents; + double b4eff = ((double) b4_count)/nevents; + double b1b3eff = b1b3_count/nevents; + double b1b4eff = b1b4_count/nevents; + double b3b4eff = b3b4_count/nevents; + double b1b3b4eff = b1b3b4_count/nevents; + TParameter break1_eff("Breakup1 Efficiency", b1eff); + TParameter break3_eff("Breakup3 Efficiency", b3eff); + TParameter break4_eff("Breakup4 Efficiency", b4eff); + TParameter b1b3_eff("Break1 Coincident with Break3", b1b3eff); + TParameter b1b4_eff("Break1 Coincident with Break4", b1b4eff); + TParameter b3b4_eff("Break3 Coincident with Break4", b3b4eff); + TParameter b1b3b4_eff("All Breakups Coincident", b1b3b4eff); + + input->cd(); + break1_eff.Write(); + break3_eff.Write(); + break4_eff.Write(); + b1b3_eff.Write(); + b1b4_eff.Write(); + b3b4_eff.Write(); + b1b3b4_eff.Write(); + table->Write(); + input->Close(); +} \ No newline at end of file diff --git a/src/AngularDistribution.cpp b/src/AngularDistribution.cpp new file mode 100644 index 0000000..453b862 --- /dev/null +++ b/src/AngularDistribution.cpp @@ -0,0 +1,104 @@ +#include "AngularDistribution.h" +#include +#include +#include +#include "LegendrePoly.h" + +AngularDistribution::AngularDistribution() : + generator(nullptr), branchingRatio(1.0), L(0), isoFlag(true) +{ +} + +AngularDistribution::AngularDistribution(const std::string& file) : + generator(nullptr), branchingRatio(1.0), L(0), isoFlag(true) +{ + ReadDistributionFile(file); +} + +AngularDistribution::~AngularDistribution() {} + +void AngularDistribution::ReadDistributionFile(const std::string& file) { + + if(file == "none" || file == "") { + L=0; + branchingRatio=1.0; + constants.clear(); + constants.push_back(0.5); + isoFlag = true; + return; + } + + std::ifstream input(file); + std::string junk; + int l; + double par; + + if(!input.is_open()) { + std::cerr<<"Unable to open distribution file. All values reset to default."<>junk>>l; + while(input>>junk) { + input>>par; + constants.push_back(par); + } + input.close(); + + if(constants.size() != ((unsigned int) l+1)) { + std::cerr<<"Unexpected number of constants for given angular momentum! Expected "<Uniform(-1.0, 1.0); + + double test, probability; + double costheta; + + test = generator->Uniform(0.0, 1.0); + if(test > branchingRatio) return -10; + + do { + probability = 0.0; + costheta = generator->Uniform(-1.0, 1.0); + test = generator->Uniform(0.0, 1.0); + for(unsigned int i=0; i probability); + + return costheta; +} \ No newline at end of file diff --git a/src/EnergyLoss.cpp b/src/EnergyLoss.cpp index 4e970c1..266d942 100644 --- a/src/EnergyLoss.cpp +++ b/src/EnergyLoss.cpp @@ -29,7 +29,7 @@ void EnergyLoss::SetTargetComponents(std::vector& Zt, std::vector& At, for(unsigned int i=0; i MAX_Z) { - throw ELossException(); + throw ELossException("Maximum allowed target Z exceeded"); } } targ_composition.resize(Stoich.size()); @@ -136,7 +136,7 @@ double EnergyLoss::GetElectronicStoppingPower(double energy) { double e_per_u = energy/MP; std::vector values; if(e_per_u > MAX_H_E_PER_U) { - throw ELossException(); + throw ELossException("Exceeded maximum allowed energy per nucleon"); } else if (e_per_u > 1000.0) { for(auto& z: ZT) { values.push_back(Hydrogen_dEdx_High(e_per_u, energy, z)); @@ -150,11 +150,11 @@ double EnergyLoss::GetElectronicStoppingPower(double energy) { values.push_back(Hydrogen_dEdx_Low(e_per_u, z)); } } else { - throw ELossException(); + throw ELossException("Negative energy per nucleon"); } if(values.size() == 0) { - throw ELossException(); + throw ELossException("Size of value array is 0. Unable to iterate over target components"); } if(ZP > 1) { //not hydrogen, need to account for effective charge @@ -192,6 +192,9 @@ double EnergyLoss::GetNuclearStoppingPower(double energy) { /*Wrapper function for aquiring total stopping (elec + nuc)*/ double EnergyLoss::GetTotalStoppingPower(double energy) { + if(ZP == 0) { + return GetNuclearStoppingPower(energy); + } return GetElectronicStoppingPower(energy)+GetNuclearStoppingPower(energy); } diff --git a/src/Kinematics.cpp b/src/Kinematics.cpp index 0126888..b51de2a 100644 --- a/src/Kinematics.cpp +++ b/src/Kinematics.cpp @@ -1,4 +1,5 @@ #include "Kinematics.h" +#include "MaskFile.h" #include #include @@ -109,25 +110,31 @@ bool Kinematics::LoadConfig(const char* filename) { input>>junk; } - double par1, par2, L1, L2; + double par1, par2; + std::string dfile1, dfile2; getline(input, junk); getline(input, junk); input>>junk>>m_nsamples; input>>junk>>par1>>junk>>par2; sys->SetBeamDistro(par1, par2); + input>>junk>>par1; + sys->SetReactionThetaType(par1); input>>junk>>par1>>junk>>par2; sys->SetTheta1Range(par1, par2); input>>junk>>par1>>junk>>par2; + sys->SetPhi1Range(par1, par2); + input>>junk>>par1>>junk>>par2; sys->SetExcitationDistro(par1, par2); - input>>junk>>L1; - input>>junk>>L2; - sys->SetDecay1AngularMomentum(L1); - sys->SetDecay2AngularMomentum(L2); + input>>junk>>dfile1; + input>>junk>>dfile2; + sys->SetDecay1Distribution(dfile1); + sys->SetDecay2Distribution(dfile2); sys->SetRandomGenerator(global_generator); std::cout<<"Reaction equation: "<GetDecay1AngularMomentum()<<" Decay2 angular momentum: "<GetDecay2AngularMomentum()<GetDecay1BranchingRatio()<<" Decay2 total branching ratio: "<GetDecay2BranchingRatio()<Branch("ejectile", &eject); tree->Branch("residual", &residual); } + //For progress tracking int percent5 = 0.05*m_nsamples; int count = 0; int npercent = 0; + std::vector data; + data.resize(4); for(int i=0; iRunSystem(); + if(save_tree_flag) { targ = ConvertNucleus(sys->GetTarget()); proj = ConvertNucleus(sys->GetProjectile()); @@ -217,8 +229,10 @@ void Kinematics::RunOneStepRxn() { plotman.FillData(sys->GetEjectile()); plotman.FillData(sys->GetResidual()); } + } + output->cd(); if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite); if(do_plotter_flag) { @@ -226,6 +240,7 @@ void Kinematics::RunOneStepRxn() { plotman.ClearTable(); } output->Close(); + } void Kinematics::RunOneStepDecay() { @@ -264,6 +279,7 @@ void Kinematics::RunOneStepDecay() { plotman.FillData(sys->GetEjectile()); plotman.FillData(sys->GetResidual()); } + } output->cd(); @@ -282,6 +298,7 @@ void Kinematics::RunTwoStep() { return; } + TFile* output = TFile::Open(m_outfile_name.c_str(), "RECREATE"); TTree* tree; NucData targ, proj, eject, residual, break1, break2; @@ -294,6 +311,8 @@ void Kinematics::RunTwoStep() { tree->Branch("breakup1", &break1); tree->Branch("breakup2", &break2); } + + //For progress tracking int percent5 = 0.05*m_nsamples; @@ -325,8 +344,10 @@ void Kinematics::RunTwoStep() { plotman.FillData(this_sys->GetBreakup1()); plotman.FillData(this_sys->GetBreakup2()); } + } + output->cd(); if(save_tree_flag) tree->Write(tree->GetName(), TObject::kOverwrite); if(do_plotter_flag) { @@ -334,6 +355,7 @@ void Kinematics::RunTwoStep() { plotman.ClearTable(); } output->Close(); + } void Kinematics::RunThreeStep() { diff --git a/src/LegendrePoly.cpp b/src/LegendrePoly.cpp index ddeeace..38d7914 100644 --- a/src/LegendrePoly.cpp +++ b/src/LegendrePoly.cpp @@ -1,4 +1,5 @@ #include "LegendrePoly.h" +#include double P_l(int l, double x) { if(l == 0) { @@ -6,6 +7,26 @@ double P_l(int l, double x) { } else if (l == 1) { return x; } else { - return ((2.0*l + 1.0)*x*P_l(l-1, x) - (l-1.0)*P_l(l-2, x))/(double(l)); + return (2.0*l - 1.0)/l*x*P_l(l-1, x) - (l-1.0)/l*P_l(l-2, x); } +} + +double Normed_P_l_sq(int l, double x) { + return (2.0*l+1.0)/2.0*std::pow(P_l(l, x), 2.0); +} + +double P_0(double x) { + return 1.0; +} + +double P_1(double x) { + return x; +} + +double P_2(double x) { + return 0.5*(3.0*x*x -1.0); +} + +double P_l_ROOT(double* x, double* pars) { + return P_l(pars[0], x[0]); } \ No newline at end of file diff --git a/src/MaskFile.cpp b/src/MaskFile.cpp new file mode 100644 index 0000000..af9b279 --- /dev/null +++ b/src/MaskFile.cpp @@ -0,0 +1,179 @@ +#include "MaskFile.h" +#include +#include +#include + +namespace Mask { + +MaskFile::MaskFile() : + file_type(MaskFile::FileType::none), filename(""), buffer_position(0), count(0), file() +{ +} + +MaskFile::MaskFile(const std::string& name, MaskFile::FileType type) : + file_type(type), filename(name), buffer_position(0), count(0), file() +{ + Open(filename, type); +} + +bool MaskFile::Open(const std::string& name, MaskFile::FileType type) { + if(IsOpen()) { + std::cerr<<"Attempted to open file that is already open!"< 0) { + FlushBuffer(); + } + file.close(); +} + +void MaskFile::WriteHeader(std::vector& nuclei) { + file<& data) { + if(count == 0 && data_buffer.size() == 0) { + WriteHeader(data); + } + + data_buffer.push_back(data); + if(data_buffer.size() == buffersize) { + FlushBuffer(); + } +} + +void MaskFile::FlushBuffer() { + for(auto& event : data_buffer) { + file< MaskFile::ReadData() { + if(buffer_position == data_buffer.size()) { + FillBuffer(); + buffer_position = 0; + } + std::vector data = data_buffer[buffer_position]; + buffer_position++; + + return data; +} + +void MaskFile::FillBuffer() { + std::vector data; + std::string line; + std::stringstream linebuf; + std::string junk; + int Z, A; + double E, p, theta, phi, tcm; + data_buffer.clear(); + + while(data_buffer.size() <= buffersize) { + if(!std::getline(file, line)) break; + linebuf.str(line); + while(linebuf>>Z) { + linebuf>>A; + linebuf>>E>>junk>>junk>>p>>theta>>phi>>tcm; + data.emplace_back(Z, A); + data[data.size()-1].SetVectorSpherical(theta, phi, p, E); + data[data.size()-1].SetThetaCM(tcm); + } + data_buffer.push_back(data); + data.clear(); + } +} + +}; \ No newline at end of file diff --git a/src/Plotter.cpp b/src/Plotter.cpp index dfe8df2..f0b233f 100644 --- a/src/Plotter.cpp +++ b/src/Plotter.cpp @@ -25,12 +25,12 @@ void Plotter::FillData(const Nucleus& nuc) { std::string ex_name = sym + "_ex"; std::string ex_title = ex_name + ";E_{ex} (MeV);counts"; std::string angdist_name = sym+"_angDist"; - std::string angdist_title = angdist_name+";#theta_{CM};counts"; + std::string angdist_title = angdist_name+";cos#right(#theta_{CM}#left);counts"; MyFill(ke_vs_th_name.c_str(), ke_vs_th_title.c_str(), nuc.GetTheta()*rad2deg, nuc.GetKE(), 2); MyFill(ke_vs_ph_name.c_str(), ke_vs_ph_title.c_str(), nuc.GetPhi()*rad2deg, nuc.GetKE(), 4); MyFill(ex_name.c_str(),ex_title.c_str(),260,-1.0,25,nuc.GetExcitationEnergy()); - MyFill(angdist_name.c_str(), angdist_title.c_str(),180,0,180,nuc.GetThetaCM()*rad2deg); + MyFill(angdist_name.c_str(), angdist_title.c_str(),100,-1.0,1.0,std::cos(nuc.GetThetaCM())); } diff --git a/src/QQQDetector.cpp b/src/QQQDetector.cpp new file mode 100644 index 0000000..40d4a51 --- /dev/null +++ b/src/QQQDetector.cpp @@ -0,0 +1,155 @@ +#include "QQQDetector.h" + +QQQDetector::QQQDetector(double R_in, double R_out, double deltaPhi, double phiCentral, double z, double x, double y) : + m_Rinner(R_in), m_Router(R_out), m_deltaPhi(deltaPhi), m_phiCentral(phiCentral), m_translation(x,y,z) +{ + m_deltaR = (m_Router - m_Rinner)/nrings; + m_deltaPhi_per_wedge = m_deltaPhi/nwedges; + m_ZRot.SetAngle(m_phiCentral); + m_ringCoords.resize(nrings); + m_wedgeCoords.resize(nwedges); + for(auto& ring : m_ringCoords) { + ring.resize(4); + } + for(auto& wedge : m_wedgeCoords) { + wedge.resize(4); + } + + CalculateCorners(); +} + +QQQDetector::~QQQDetector() {} + +void QQQDetector::CalculateCorners() { + double x0, x1, x2, x3; + double y0, y1, y2, y3; + double z0, z1, z2, z3; + + //Generate flat ring corner coordinates + for(int i=0; i min_rho) { + for(auto& wedge : m_wedgeCoords) { + min_phi = wedge[0].GetPhi(); + max_phi = wedge[3].GetPhi(); + if(phi < min_phi && phi < max_phi) { + result.SetVectorSpherical(r_traj, theta, phi); + break; + } + } + } + } + + return result; + +} + +std::pair QQQDetector::GetTrajectoryRingWedge(double theta, double phi) { + double z_to_detector = m_translation.GetZ(); + double rho_traj = z_to_detector*std::tan(theta); + double min_rho, max_rho, min_phi, max_phi; + + + for(int r=0; r min_rho) { + for(int w=0; w min_phi && phi < max_phi) { + return std::make_pair(r, w); + } + } + } + } + + return std::make_pair(-1, -1); + +} + +Mask::Vec3 QQQDetector::GetHitCoordinates(int ringch, int wedgech) { + if(!CheckChannel(ringch) || !CheckChannel(wedgech)) { + return Mask::Vec3(); + } + + double r_center = m_Rinner + (0.5+ringch)*m_deltaR; + double phi_center = -m_deltaPhi/2.0 + (0.5+wedgech)*m_deltaPhi_per_wedge; + double x = r_center*std::cos(phi_center); + double y = r_center*std::sin(phi_center); + double z = 0; + + Mask::Vec3 hit(x, y, z); + + return TransformCoordinates(hit); +} \ No newline at end of file diff --git a/src/Reaction.cpp b/src/Reaction.cpp index 92621a0..3aac0cf 100644 --- a/src/Reaction.cpp +++ b/src/Reaction.cpp @@ -12,12 +12,12 @@ namespace Mask { Reaction::Reaction() : - target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), nuc_initFlag(false), resid_elossFlag(false) + target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), nuc_initFlag(false), resid_elossFlag(false) { } Reaction::Reaction(int zt, int at, int zp, int ap, int ze, int ae) : - target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), resid_elossFlag(false) + target(nullptr), m_bke(0), m_theta(0), m_phi(0), m_ex(0), rxnLayer(0), m_eject_theta_type(lab), resid_elossFlag(false) { SetNuclei(zt, at, zp, ap, ze, ae); } @@ -77,16 +77,28 @@ void Reaction::SetBeamKE(double bke) { m_bke = bke - target->GetProjectileEnergyLoss(reactants[1].GetZ(), reactants[1].GetA(), bke, rxnLayer, 0); }; +void Reaction::SetEjectileThetaType(int type) { + if(decayFlag) return; + if(type != center_of_mass && type != lab) return; + + m_eject_theta_type = type; +} + //Methods given by Iliadis in Nuclear Physics of Stars, Appendix C -void Reaction::CalculateReaction() { - //Target assumed at rest, with 0 excitation energy +//For use with lab frame restricted angles. May not give appropriate disribution for ejectile +void Reaction::CalculateReactionThetaLab() { reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass()); double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass())); double beam_E = m_bke + reactants[1].GetGroundStateMass(); reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E); - double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex); + + double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()); + if(m_bke < Ethresh) { + throw EnergyThresholdException(); + } + double term1 = sqrt(reactants[1].GetGroundStateMass()*reactants[2].GetGroundStateMass()*m_bke)/(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())*cos(m_theta); double term2 = (m_bke*(reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()) + reactants[3].GetGroundStateMass()*Q)/(reactants[3].GetGroundStateMass() + reactants[2].GetGroundStateMass()); double sqrt_pos_ejectKE = term1 + sqrt(term1*term1 + term2); @@ -97,7 +109,6 @@ void Reaction::CalculateReaction() { } else { ejectKE = sqrt_neg_ejectKE*sqrt_neg_ejectKE; } - double ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass())); double ejectE = ejectKE + reactants[2].GetGroundStateMass(); @@ -105,12 +116,59 @@ void Reaction::CalculateReaction() { reactants[3] = reactants[0] + reactants[1] - reactants[2]; - //energy loss for ejectile (after reaction!) ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, m_theta); ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass())); ejectE = ejectKE + reactants[2].GetGroundStateMass(); reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP, ejectE); + if(resid_elossFlag) { + double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta()); + double residP = std::sqrt(residKE*(residKE + 2.0*reactants[3].GetInvMass())); + double residE = residKE + reactants[3].GetInvMass(); + reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE); + } +} + +//Methods from original ANASEN. Gives proper distribution for inverse kinematics. +void Reaction::CalculateReactionThetaCM() { + //Target assumed at rest, with 0 excitation energy + reactants[0].SetVectorCartesian(0.,0.,0.,reactants[0].GetGroundStateMass()); + double beam_pz = std::sqrt(m_bke*(m_bke + 2.0 * reactants[1].GetGroundStateMass())); + double beam_E = m_bke + reactants[1].GetGroundStateMass(); + reactants[1].SetVectorCartesian(0.,0.,beam_pz,beam_E); + + + double Q = reactants[0].GetGroundStateMass() + reactants[1].GetGroundStateMass() - (reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() + m_ex); + + double Ethresh = -Q*(reactants[2].GetGroundStateMass()+reactants[3].GetGroundStateMass())/(reactants[2].GetGroundStateMass() + reactants[3].GetGroundStateMass() - reactants[1].GetGroundStateMass()); + if(m_bke < Ethresh) { + throw EnergyThresholdException(); + } + + auto parent = reactants[0] + reactants[1]; + double boost2lab[3]; + double boost2cm[3]; + const double* boost = parent.GetBoost(); + for(int i=0; i<3; i++) { + boost2lab[i] = boost[i]; + boost2cm[i] = boost[i]*(-1.0); + } + parent.ApplyBoost(boost2cm); + double ejectE_cm = (std::pow(reactants[2].GetGroundStateMass(), 2.0) - std::pow(reactants[3].GetGroundStateMass() + m_ex, 2.0) + std::pow(parent.GetE(),2.0))/(2.0*parent.GetE()); + double ejectP_cm = std::sqrt(ejectE_cm*ejectE_cm - std::pow(reactants[2].GetGroundStateMass(), 2.0)); + reactants[2].SetVectorSpherical(m_theta, m_phi, ejectP_cm, ejectE_cm); + reactants[2].ApplyBoost(boost2lab); + reactants[3] = reactants[0] + reactants[1] - reactants[2]; + + double ejectKE = reactants[2].GetKE(); + double ejectP = reactants[2].GetP(); + double ejectE = reactants[2].GetE(); + //energy loss for ejectile (after reaction!) + ejectKE -= target->GetEjectileEnergyLoss(reactants[2].GetZ(), reactants[2].GetA(), ejectKE, rxnLayer, reactants[2].GetTheta()); + ejectP = std::sqrt(ejectKE*(ejectKE + 2.0 * reactants[2].GetGroundStateMass())); + ejectE = ejectKE + reactants[2].GetGroundStateMass(); + reactants[2].SetVectorSpherical(reactants[2].GetTheta(), reactants[2].GetPhi(), ejectP, ejectE); + //if on, get eloss for residual (after reaction!) if(resid_elossFlag) { double residKE = reactants[3].GetKE() - target->GetEjectileEnergyLoss(reactants[3].GetZ(), reactants[3].GetA(), reactants[3].GetKE(), rxnLayer, reactants[3].GetTheta()); @@ -118,7 +176,21 @@ void Reaction::CalculateReaction() { double residE = residKE + reactants[3].GetInvMass(); reactants[3].SetVectorSpherical(reactants[3].GetTheta(), reactants[3].GetPhi(), residP, residE); } +} +void Reaction::CalculateReaction() { + switch(m_eject_theta_type) { + case center_of_mass: + { + CalculateReactionThetaCM(); + break; + } + case lab: + { + CalculateReactionThetaLab(); + break; + } + } } //Calculate in CM, where decay is isotropic diff --git a/src/ReactionSystem.cpp b/src/ReactionSystem.cpp index dfc2f7e..3cf53d1 100644 --- a/src/ReactionSystem.cpp +++ b/src/ReactionSystem.cpp @@ -5,12 +5,12 @@ namespace Mask { ReactionSystem::ReactionSystem() : - m_beamDist(0,0), m_theta1Range(0,0), m_exDist(0,0), generator(nullptr), target_set_flag(false), gen_set_flag(false), rxnLayer(0), m_sys_equation("") + m_beamDist(0,0), m_theta1Range(0,0), m_phi1Range(0,0), m_exDist(0,0), generator(nullptr), target_set_flag(false), gen_set_flag(false), rxnLayer(0), m_sys_equation("") { } ReactionSystem::ReactionSystem(std::vector& z, std::vector& a) : - m_beamDist(0,0), m_theta1Range(0,0), m_exDist(0,0), generator(nullptr), target_set_flag(false), gen_set_flag(false), rxnLayer(0), m_sys_equation("") + m_beamDist(0,0), m_theta1Range(0,0), m_phi1Range(0,0), m_exDist(0,0), generator(nullptr), target_set_flag(false), gen_set_flag(false), rxnLayer(0), m_sys_equation("") { SetNuclei(z, a); } @@ -29,6 +29,13 @@ bool ReactionSystem::SetNuclei(std::vector&z, std::vector& a) { return true; } +void ReactionSystem::SetRandomGenerator(TRandom3* gen) { + generator = gen; + decay1dist.AttachRandomNumberGenerator(gen); + decay2dist.AttachRandomNumberGenerator(gen); + gen_set_flag = true; +} + void ReactionSystem::AddTargetLayer(std::vector& zt, std::vector& at, std::vector& stoich, double thickness) { target.AddLayer(zt, at, stoich, thickness); } @@ -55,22 +62,6 @@ void ReactionSystem::SetSystemEquation() { m_sys_equation += step1.GetResidual().GetIsotopicSymbol(); } -double ReactionSystem::GetDecayTheta(int L) { - if(!gen_set_flag) return 0.0; - - double probability = 0.0; - double costheta = 0.0; - double check = 0.0; - - do { - costheta = generator->Uniform(-1.0, 1.0); - check = generator->Uniform(0.0, 1.0); - probability = std::pow(P_l(L, costheta), 2.0); - } while(check > probability); - - return std::acos(costheta); -} - void ReactionSystem::RunSystem() { if(!gen_set_flag) return; @@ -80,7 +71,13 @@ void ReactionSystem::RunSystem() { } if(step1.IsDecay()) { - double rxnTheta = GetDecayTheta(L1); + double decaycostheta = decay1dist.GetRandomCosTheta(); + if(decaycostheta == -10.0) { + step1.ResetEjectile(); + step1.ResetResidual(); + return; + } + double rxnTheta = std::acos(decay1dist.GetRandomCosTheta()); double rxnPhi = generator->Uniform(0, 2.0*M_PI); step1.SetPolarRxnAngle(rxnTheta); step1.SetAzimRxnAngle(rxnPhi); @@ -91,7 +88,7 @@ void ReactionSystem::RunSystem() { //Sample parameters double bke = generator->Gaus(m_beamDist.first, m_beamDist.second); double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second))); - double rxnPhi = 0; + double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second); double residEx = generator->Gaus(m_exDist.first, m_exDist.second); step1.SetBeamKE(bke); diff --git a/src/SabreEfficiency.cpp b/src/SabreEfficiency.cpp index 662b298..88aa55b 100644 --- a/src/SabreEfficiency.cpp +++ b/src/SabreEfficiency.cpp @@ -1,5 +1,4 @@ #include "SabreEfficiency.h" -#include "Kinematics.h" #include #include #include @@ -10,9 +9,10 @@ #include #include #include +#include SabreEfficiency::SabreEfficiency() : - m_rxn_type(-1), deadlayer(DEADLAYER_THIN) + m_rxn_type(-1), deadlayer(DEADLAYER_THIN), sabre_eloss(SABRE_THICKNESS) { detectors.reserve(5); detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI0*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); @@ -21,59 +21,39 @@ SabreEfficiency::SabreEfficiency() : detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); - Mask::Vec3 coords; - for(int i=0; i<5; i++) { - for(int j=0; j dead_z = {14}; std::vector dead_a = {28}; std::vector dead_stoich = {1}; deadlayer.SetElements(dead_z, dead_a, dead_stoich); + sabre_eloss.SetElements(dead_z, dead_a, dead_stoich); } SabreEfficiency::~SabreEfficiency() {} -void SabreEfficiency::MyFill(THashTable* table, const char* name, const char* title, int bins, float min, float max, double val) { - TH1F* h = (TH1F*) table->FindObject(name); +void SabreEfficiency::MyFill(THashTable* table, const std::string& name, const std::string& title, int bins, float min, float max, double val) { + TH1F* h = (TH1F*) table->FindObject(name.c_str()); if(h) { h->Fill(val); } else { - h = new TH1F(name, title, bins, min, max); + h = new TH1F(name.c_str(), title.c_str(), bins, min, max); h->Fill(val); table->Add(h); } } -void SabreEfficiency::MyFill(THashTable* table, const char* name, const char* title, int binsx, float minx, float maxx, int binsy, float miny, float maxy, double valx, double valy) { - TH2F* h = (TH2F*) table->FindObject(name); +void SabreEfficiency::MyFill(THashTable* table, const std::string& name, const std::string& title, int binsx, float minx, float maxx, double valx, int binsy, float miny, float maxy, double valy) { + TH2F* h = (TH2F*) table->FindObject(name.c_str()); if(h) { h->Fill(valx, valy); } else { - h = new TH2F(name, title, binsx, minx, maxx, binsy, miny, maxy); + h = new TH2F(name.c_str(), title.c_str(), binsx, minx, maxx, binsy, miny, maxy); h->Fill(valx, valy); table->Add(h); } } -void SabreEfficiency::CalculateEfficiency(const char* file) { +void SabreEfficiency::CalculateEfficiency(const std::string& file) { std::cout<<"----------SABRE Efficiency Calculation----------"<Get("DataTree"); - - Mask::NucData* eject = new Mask::NucData(); - Mask::NucData* resid = new Mask::NucData(); - - - tree->SetBranchAddress("ejectile", &eject); - tree->SetBranchAddress("residual", &resid); - - double nevents = tree->GetEntries(); - std::vector resid_xs, eject_xs; - std::vector resid_ys, eject_ys; - std::vector resid_zs, eject_zs; - - //Progress tracking - int percent5 = nevents*0.05; - int count = 0; - int npercent = 0; +void SabreEfficiency::DrawDetectorSystem(const std::string& filename) { + TFile* output = TFile::Open(filename.c_str(), "RECREATE"); + std::vector ringxs, ringys, ringzs; + std::vector wedgexs, wedgeys, wedgezs; Mask::Vec3 coords; - double thetaIncident, eloss; + for(int i=0; i<5; i++) { + for(int j=0; jGetEntries(); i++) { - if(++count == percent5) {//Show progress every 5% - npercent++; - count = 0; - std::cout<<"\rPercent completed: "<GetEntry(i); - - if(eject->KE >= ENERGY_THRESHOLD) { - for(int j=0; j<5; j++) { - auto& det = detectors[j]; - auto chan = det.GetTrajectoryRingWedge(eject->theta, eject->phi); - if(chan.first != -1 && chan.second != -1) { - if(dmap.IsDead(j, chan.first, 0) || dmap.IsDead(j, chan.second, 1)) break; //dead channel check - coords = det.GetTrajectoryCoordinates(eject->theta, eject->phi); - thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR())); - eloss = deadlayer.getEnergyLossTotal(eject->Z, eject->A, eject->KE, M_PI - thetaIncident); - if((eject->KE - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check - - eject_xs.push_back(coords.GetX()); - eject_ys.push_back(coords.GetY()); - eject_zs.push_back(coords.GetZ()); - break; - } - } - } - - if(resid->KE > ENERGY_THRESHOLD) { - for(int j=0; j<5; j++) { - auto& det = detectors[j]; - auto chan = det.GetTrajectoryRingWedge(resid->theta, resid->phi); - if(chan.first != -1 && chan.second != -1) { - if(dmap.IsDead(j, chan.first, 0) || dmap.IsDead(j, chan.second, 1)) break; - coords = det.GetTrajectoryCoordinates(resid->theta, resid->phi); - thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR())); - eloss = deadlayer.getEnergyLossTotal(resid->Z, resid->A, resid->KE, M_PI - thetaIncident); - if((resid->KE - eloss) <= ENERGY_THRESHOLD) break; - - resid_xs.push_back(coords.GetX()); - resid_ys.push_back(coords.GetY()); - resid_zs.push_back(coords.GetZ()); - break; - } - } - } - - } - - double ejecteff = ((double) eject_xs.size())/nevents; - double resideff = ((double) resid_xs.size())/nevents; - TParameter eject_eff("Light Breakup Efficiency", ejecteff); - TParameter resid_eff("Heavy Breakup Efficiency", resideff); - - TGraph2D* gde = new TGraph2D(eject_xs.size(), &(eject_xs[0]), &(eject_ys[0]), &(eject_zs[0])); - gde->SetName("detected_eject_points"); - gde->SetMarkerStyle(1); - gde->SetMarkerColor(2); - - TGraph2D* gr = new TGraph2D(ringxs.size(), &(ringxs[0]), &(ringys[0]), &(ringzs[0])); + TGraph2D* gr = new TGraph2D(ringxs.size(), &(ringxs[0]), &(ringys[0]), &(ringzs[0])); gr->SetName("ring_detector_edges"); gr->SetTitle("SABRE Detector; x(m); y(m); z(m)"); gr->SetMarkerStyle(1); @@ -201,26 +126,127 @@ void SabreEfficiency::RunDecay(const char* file) { gw->SetMarkerStyle(1); TCanvas* canvas = new TCanvas(); - canvas->SetName("detectors_and_particles"); + canvas->SetName("detector_system"); canvas->cd(); gr->Draw("AP"); gw->Draw("same P"); - gde->Draw("same P"); + + canvas->Write(); + gr->Write(); + gw->Write(); + + output->Close(); +} + +double SabreEfficiency::RunConsistencyCheck() { + double theta, phi; + double npoints = 5.0*16.0*4.0; + int count=0; + for(int h=0; h<5; h++) { + for(int j=0; j<16; j++) { + for(int k=0; k<4; k ++) { + theta = detectors[h].GetRingTiltCoords(j, k).GetTheta(); + phi = detectors[h].GetRingTiltCoords(j, k).GetPhi(); + for(int i=0; i<5; i++) { + auto channels = detectors[i].GetTrajectoryRingWedge(theta, phi); + if(channels.first != -1) { + count++; + } + } + } + } + } + + return ((double)count)/npoints; +} + +/*Returns if detected, as well as total energy deposited in SABRE*/ +std::pair SabreEfficiency::IsSabre(Mask::NucData* nucleus) { + if(nucleus->KE <= ENERGY_THRESHOLD) { + return std::make_pair(false, 0.0); + } + + Mask::Vec3 coords; + double thetaIncident, eloss, e_deposited; + for(int i=0; i<5; i++) { + auto chan = detectors[i].GetTrajectoryRingWedge(nucleus->theta, nucleus->phi); + if(chan.first != -1 && chan.second != -1) { + if(dmap.IsDead(i, chan.first, 0) || dmap.IsDead(i, chan.second, 1)) break; //dead channel check + coords = detectors[i].GetTrajectoryCoordinates(nucleus->theta, nucleus->phi); + thetaIncident = std::acos(coords.Dot(detectors[i].GetNormTilted())/(coords.GetR())); + eloss = deadlayer.getEnergyLossTotal(nucleus->Z, nucleus->A, nucleus->KE, M_PI - thetaIncident); + if((nucleus->KE - eloss) <= ENERGY_THRESHOLD) break; //deadlayer check + e_deposited = sabre_eloss.getEnergyLossTotal(nucleus->Z, nucleus->A, nucleus->KE - eloss, M_PI - thetaIncident); + return std::make_pair(true, e_deposited); + } + } + + return std::make_pair(false,0.0); +} + +void SabreEfficiency::RunDecay(const std::string& filename) { + TFile* input = TFile::Open(filename.c_str(), "UPDATE"); + TTree* tree = (TTree*) input->Get("DataTree"); + + THashTable* table = new THashTable(); + + Mask::NucData* eject = new Mask::NucData(); + Mask::NucData* resid = new Mask::NucData(); + + + tree->SetBranchAddress("ejectile", &eject); + tree->SetBranchAddress("residual", &resid); + + double nevents = tree->GetEntries(); + + //Progress tracking + int percent5 = nevents*0.05; + int count = 0; + int npercent = 0; + + std::pair eject_result, resid_result; + int detected_eject=0, detected_resid=0; + for(int i=0; iGetEntries(); i++) { + if(++count == percent5) {//Show progress every 5% + npercent++; + count = 0; + std::cout<<"\rPercent completed: "<GetEntry(i); + + eject_result = IsSabre(eject); + resid_result = IsSabre(resid); + + if(eject_result.first) { + detected_eject++; + MyFill(table, "detected_eject_theta_ke","detected ejectile theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,eject->theta/DEG2RAD,300,0.0,30.0,eject->KE); + MyFill(table, "detected_eject_theta_edep","detected ejectile theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,eject->theta/DEG2RAD,300,0.0,30.0,eject_result.second); + } + + if(resid_result.first) { + detected_resid++; + MyFill(table, "detected_resid_theta_ke","detected residual theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,resid->theta/DEG2RAD,300,0.0,30.0,resid->KE); + MyFill(table, "detected_resid_theta_edep","detected residual theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,resid->theta/DEG2RAD,300,0.0,30.0,resid_result.second); + } + } + + double ejecteff = ((double) detected_eject)/nevents; + double resideff = ((double) detected_resid)/nevents; + TParameter eject_eff("Ejectile Efficiency", ejecteff); + TParameter resid_eff("Residual Efficiency", resideff); input->cd(); eject_eff.Write(); resid_eff.Write(); - gr->Write(); - gw->Write(); - gde->Write(); - canvas->Write(); + table->Write(); input->Close(); } -void SabreEfficiency::Run2Step(const char* file) { +void SabreEfficiency::Run2Step(const std::string& filename) { - TFile* input = TFile::Open(file, "UPDATE"); + TFile* input = TFile::Open(filename.c_str(), "UPDATE"); TTree* tree = (TTree*) input->Get("DataTree"); Mask::NucData* break1 = new Mask::NucData(); @@ -232,20 +258,17 @@ void SabreEfficiency::Run2Step(const char* file) { tree->SetBranchAddress("breakup2", &break2); double nevents = tree->GetEntries(); - std::vector b1_thetas, b2_thetas; - std::vector b1_phis, b2_phis; - std::vector b1_kes, b2_kes; //Progress tracking int percent5 = nevents*0.05; int count = 0; int npercent = 0; - int cnt_09=0, cnt_08=0, cnt_07=0, cnt_06=0; double costheta_cm =0; - Mask::Vec3 coords; - double thetaIncident, eloss; + std::pair break1_result, break2_result; + int detected_b1=0, detected_b2=0; + for(int i=0; iGetEntries(); i++) { if(++count == percent5) {//Show progress every 5% npercent++; @@ -255,97 +278,44 @@ void SabreEfficiency::Run2Step(const char* file) { tree->GetEntry(i); - if(break1->KE >= ENERGY_THRESHOLD) { - for(int j=0; j<5; j++) { - auto& det = detectors[j]; - auto chan = det.GetTrajectoryRingWedge(break1->theta, break1->phi); - if(chan.first != -1 && chan.second != -1) { - if(dmap.IsDead(j, chan.first, 0) || dmap.IsDead(j, chan.second, 1)) break; - costheta_cm = std::cos(break1->theta_cm); - double this_bin; - std::string this_name; - for(int k=-9; k<=10; k++) { - this_bin = 0.1*k; - this_name = "counter_histo_bin"+std::to_string(k); - if(costheta_cm > this_bin-0.05 && costheta_cm < this_bin+0.05) { - MyFill(table, this_name.c_str(), "ctheta",20,-1.0,1.0,costheta_cm); - break; - } - } - /* - if(costheta_cm > -0.95 && costheta_cm < -0.85) { - cnt_09++; - } else if(costheta_cm > -0.85 && costheta_cm < -0.75) { - cnt_08++; - } else if(costheta_cm > -0.75 && costheta_cm < -0.65) { - cnt_07++; - } else if(costheta_cm > -0.65 && costheta_cm < -0.55) { - cnt_06++; - } - */ - coords = det.GetTrajectoryCoordinates(break1->theta, break1->phi); - thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR())); - eloss = deadlayer.getEnergyLossTotal(break1->Z, break1->A, break1->KE, M_PI - thetaIncident); - if((break1->KE - eloss) <= ENERGY_THRESHOLD) break; - + break1_result = IsSabre(break1); + break2_result = IsSabre(break2); - b1_thetas.push_back(break1->theta); - b1_phis.push_back(break1->phi); - b1_kes.push_back(break1->KE); - break; - } - } + if(break1_result.first) { + detected_b1++; + costheta_cm = std::cos(break1->theta_cm); + MyFill(table,"detected_break1_cm_theta","cos(#theta_{CM})",20,-1.0,1.0,costheta_cm); + MyFill(table, "detected_break1_theta_ke","detected break1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1->KE); + MyFill(table, "detected_break1_theta_edep","detected break1 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1_result.second); } - if(break2->KE > ENERGY_THRESHOLD) { - for(int j=0; j<5; j++) { - auto& det = detectors[j]; - auto chan = det.GetTrajectoryRingWedge(break2->theta, break2->phi); - if(chan.first != -1 && chan.second != -1) { - if(dmap.IsDead(j, chan.first, 0) || dmap.IsDead(j, chan.second, 1)) break; - coords = det.GetTrajectoryCoordinates(break2->theta, break2->phi); - thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR())); - eloss = deadlayer.getEnergyLossTotal(break2->Z, break2->A, break2->KE, M_PI - thetaIncident); - if((break2->KE - eloss) <= ENERGY_THRESHOLD) break; - - b2_thetas.push_back(break2->theta); - b2_phis.push_back(break2->phi); - b2_kes.push_back(break2->KE); - break; - } - } + if(break2_result.first) { + detected_b2++; + MyFill(table, "detected_break2_theta_ke","detected break2 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break2->theta/DEG2RAD,300,0.0,30.0,break2->KE); + MyFill(table, "detected_break2_theta_edep","detected break2 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break2->theta/DEG2RAD,300,0.0,30.0,break2_result.second); + } } - double b1eff = ((double) b1_thetas.size())/nevents; - double b2eff = ((double) b2_thetas.size())/nevents; - double b1eff_09 = cnt_09/nevents; - double b1eff_08 = cnt_08/nevents; - double b1eff_07 = cnt_07/nevents; - double b1eff_06 = cnt_06/nevents; - TParameter break1_eff("Light Breakup Efficiency", b1eff); - TParameter break1_eff_09("Light Breakup Efficiency CosTheta -0.9", b1eff_09); - TParameter break1_eff_08("Light Breakup Efficiency CosTheta -0.8", b1eff_08); - TParameter break1_eff_07("Light Breakup Efficiency CosTheta -0.7", b1eff_07); - TParameter break1_eff_06("Light Breakup Efficiency CosTheta -0.6", b1eff_06); - TParameter break2_eff("Heavy Breakup Efficiency", b2eff); + double b1eff = ((double) detected_b1)/nevents; + double b2eff = ((double) detected_b2)/nevents; + TParameter break1_eff("Breakup1 Efficiency", b1eff); + TParameter break2_eff("Breakup2 Efficiency", b2eff); input->cd(); table->Write(); break1_eff.Write(); - break1_eff_09.Write(); - break1_eff_08.Write(); - break1_eff_07.Write(); - break1_eff_06.Write(); break2_eff.Write(); input->Close(); } -void SabreEfficiency::Run3Step(const char* file) { - TFile* input = TFile::Open(file, "UPDATE"); +void SabreEfficiency::Run3Step(const std::string& filename) { + TFile* input = TFile::Open(filename.c_str(), "UPDATE"); TTree* tree = (TTree*) input->Get("DataTree"); + THashTable* table = new THashTable(); + Mask::NucData* break1 = new Mask::NucData(); Mask::NucData* break3 = new Mask::NucData(); Mask::NucData* break4 = new Mask::NucData(); @@ -356,23 +326,16 @@ void SabreEfficiency::Run3Step(const char* file) { tree->SetBranchAddress("breakup4", &break4); double nevents = tree->GetEntries(); - std::vector b1_thetas, b3_thetas, b4_thetas; - std::vector b1_phis, b3_phis, b4_phis; - std::vector b1_kes, b3_kes, b4_kes; //Progress tracking int percent5 = nevents*0.05; int count = 0; int npercent = 0; - bool break1_flag, break3_flag, break4_flag; + std::pair break1_result, break3_result, break4_result; + int detected_b1=0, detected_b3=0, detected_b4=0; int b1b3_count=0, b1b4_count=0, b3b4_count=0, b1b3b4_count=0; - TH2F* b3b4_thetas = new TH2F("b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0,180,180,0,180); - TH2F* b3b4_kes = new TH2F("b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0,15,150,0,15); - TH2F* b3b4_phis = new TH2F("b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0,360,360,0,360); - Mask::Vec3 coords; - double thetaIncident, eloss; for(int i=0; iGetEntries(); i++) { if(++count == percent5) {//Show progress every 5% @@ -381,100 +344,61 @@ void SabreEfficiency::Run3Step(const char* file) { std::cout<<"\rPercent completed: "<GetEntry(i); - if(break1->KE > ENERGY_THRESHOLD) { - for(int j=0; j<5; j++) { - auto& det = detectors[j]; - auto chan = det.GetTrajectoryRingWedge(break1->theta, break1->phi); - if(chan.first != -1 && chan.second != -1) { - coords = det.GetTrajectoryCoordinates(break1->theta, break1->phi); - thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR())); - eloss = deadlayer.getEnergyLossTotal(break1->Z, break1->A, break1->KE, M_PI - thetaIncident); - if((break1->KE - eloss) <= ENERGY_THRESHOLD) break; - - b1_thetas.push_back(break1->theta); - b1_phis.push_back(break1->phi); - b1_kes.push_back(break1->KE); - break1_flag = true; - break; - } - } + if(break1_result.first) { + detected_b1++; + MyFill(table, "detected_break1_theta_ke","detected break1 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1->KE); + MyFill(table, "detected_break1_theta_edep","detected break1 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break1->theta/DEG2RAD,300,0.0,30.0,break1_result.second); } - - if(break3->KE > ENERGY_THRESHOLD) { - for(int j=0; j<5; j++) { - auto& det = detectors[j]; - auto chan = det.GetTrajectoryRingWedge(break3->theta, break3->phi); - if(chan.first != -1 && chan.second != -1) { - coords = det.GetTrajectoryCoordinates(break3->theta, break3->phi); - thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR())); - eloss = deadlayer.getEnergyLossTotal(break3->Z, break3->A, break3->KE, M_PI - thetaIncident); - if((break3->KE - eloss) <= ENERGY_THRESHOLD) break; - - b3_thetas.push_back(break3->theta); - b3_phis.push_back(break3->phi); - b3_kes.push_back(break3->KE); - break3_flag = true; - break; - } - } + if(break3_result.first) { + detected_b3++; + MyFill(table, "detected_break3_theta_ke","detected break3 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break3->theta/DEG2RAD,300,0.0,30.0,break3->KE); + MyFill(table, "detected_break3_theta_edep","detected break3 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break3->theta/DEG2RAD,300,0.0,30.0,break3_result.second); } - if(break4->KE > ENERGY_THRESHOLD) { - for(int j=0; j<5; j++) { - auto& det = detectors[j]; - auto chan = det.GetTrajectoryRingWedge(break4->theta, break4->phi); - if(chan.first != -1 && chan.second != -1) { - coords = det.GetTrajectoryCoordinates(break4->theta, break4->phi); - thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR())); - eloss = deadlayer.getEnergyLossTotal(break4->Z, break4->A, break4->KE, M_PI - thetaIncident); - if((break4->KE - eloss) <= ENERGY_THRESHOLD) break; - - b4_thetas.push_back(break4->theta); - b4_phis.push_back(break4->phi); - b4_kes.push_back(break4->KE); - break4_flag = true; - break; - } - } + if(break4_result.first) { + detected_b1++; + MyFill(table, "detected_break4_theta_ke","detected break4 theta vs ke;#theta (deg);KE(MeV)",180,0.0,180.0,break4->theta/DEG2RAD,300,0.0,30.0,break4->KE); + MyFill(table, "detected_break4_theta_edep","detected break4 theta vs edep;#theta (deg);E deposited(MeV)",180,0.0,180.0,break4->theta/DEG2RAD,300,0.0,30.0,break4_result.second); } - if(break1_flag && break3_flag && break4_flag) { + if(break1_result.first && break3_result.first && break4_result.first) { b1b3b4_count++; b1b3_count++; b1b4_count++; b3b4_count++; - b3b4_thetas->Fill(b3_thetas[b3_thetas.size()-1]/DEG2RAD, b4_thetas[b4_thetas.size()-1]/DEG2RAD); - b3b4_kes->Fill(b3_kes[b3_kes.size()-1], b4_kes[b4_kes.size()-1]); - b3b4_phis->Fill(b3_phis[b3_phis.size()-1]/DEG2RAD, b4_phis[b4_phis.size()-1]/DEG2RAD); - } else if(break1_flag && break3_flag) { + MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/DEG2RAD,180,0,180,break4->theta/DEG2RAD); + MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0.0,15.0,break3->KE,150,0.0,15.0,break4->KE); + MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/DEG2RAD,360,0.0,360.0,break4->phi/DEG2RAD); + } else if(break1_result.first && break3_result.first) { b1b3_count++; - } else if(break1_flag && break4_flag) { + } else if(break1_result.first && break4_result.first) { b1b4_count++; - } else if(break3_flag && break4_flag) { + } else if(break3_result.first && break4_result.first) { b3b4_count++; - b3b4_thetas->Fill(b3_thetas[b3_thetas.size()-1]/DEG2RAD, b4_thetas[b4_thetas.size()-1]/DEG2RAD); - b3b4_kes->Fill(b3_kes[b3_kes.size()-1], b4_kes[b4_kes.size()-1]); - b3b4_phis->Fill(b3_phis[b3_phis.size()-1]/DEG2RAD, b4_phis[b4_phis.size()-1]/DEG2RAD); + MyFill(table,"b3b4_theta_theta","b3b4_theta_theta;#theta_{3};#theta_{4}",180,0.0,180.0,break3->theta/DEG2RAD,180,0,180,break4->theta/DEG2RAD); + MyFill(table,"b3b4_ke_ke","b3b4_ke_ke;KE_{3};KE_{4}",150,0.0,15.0,break3->KE,150,0.0,15.0,break4->KE); + MyFill(table,"b3b4_phi_phi","b3b4_phi_phi;#phi_{3};#phi_{4}",360,0.0,360.0,break3->phi/DEG2RAD,360,0.0,360.0,break4->phi/DEG2RAD); } } - double b1eff = ((double) b1_thetas.size())/nevents; - double b3eff = ((double) b3_thetas.size())/nevents; - double b4eff = ((double) b4_thetas.size())/nevents; + double b1eff = ((double) detected_b1)/nevents; + double b3eff = ((double) detected_b3)/nevents; + double b4eff = ((double) detected_b4)/nevents; double b1b3eff = b1b3_count/nevents; double b1b4eff = b1b4_count/nevents; double b3b4eff = b3b4_count/nevents; double b1b3b4eff = b1b3b4_count/nevents; - TParameter break1_eff("Light Initial Breakup Efficiency", b1eff); - TParameter break3_eff("Light Final Breakup Efficiency", b3eff); - TParameter break4_eff("Heavy Final Breakup Efficiency", b4eff); + TParameter break1_eff("Breakup1 Efficiency", b1eff); + TParameter break3_eff("Breakup3 Efficiency", b3eff); + TParameter break4_eff("Breakup4 Efficiency", b4eff); TParameter b1b3_eff("Break1 Coincident with Break3", b1b3eff); TParameter b1b4_eff("Break1 Coincident with Break4", b1b4eff); TParameter b3b4_eff("Break3 Coincident with Break4", b3b4eff); @@ -488,8 +412,7 @@ void SabreEfficiency::Run3Step(const char* file) { b1b4_eff.Write(); b3b4_eff.Write(); b1b3b4_eff.Write(); - b3b4_thetas->Write(); - b3b4_phis->Write(); - b3b4_kes->Write(); + table->Write(); + input->Close(); } diff --git a/src/StripDetector.cpp b/src/StripDetector.cpp new file mode 100644 index 0000000..80ad21b --- /dev/null +++ b/src/StripDetector.cpp @@ -0,0 +1,141 @@ +#include "StripDetector.h" + +/* + Corner layout for each strip in the un-rotated frame + 0--------------------------1 + | | y + | | ^ + | | | + | | | + | | | + 2--------------------------3 z<--------X + x +*/ + +StripDetector::StripDetector(int ns, double len, double wid, double cphi, double cz, double crho) : + m_random(nullptr) +{ + + num_strips = ns; + + length = std::fabs(len); + total_width = std::fabs(wid); + front_strip_width = total_width/num_strips; + back_strip_length = length/num_strips; + + while (cphi < 0) cphi += 2*M_PI; + center_phi = cphi; + center_z = cz; + center_rho = std::fabs(crho); + + zRot.SetAngle(center_phi); + + front_strip_coords.resize(num_strips); + back_strip_coords.resize(num_strips); + rotated_front_strip_coords.resize(num_strips); + rotated_back_strip_coords.resize(num_strips); + for(int i=0; iUniform(0.0, 1.0))*front_strip_width; + } else { + y = -total_width/2.0 + (front_stripch+0.5)*front_strip_width; + } + + //recall we're still assuming phi=0 det: + Mask::Vec3 coords(center_rho, y, front_strip_ratio*(length/2) + center_z); + + //NOW rotate by appropriate phi + return zRot*coords; + +} + +std::pair StripDetector::GetChannelRatio(double theta, double phi) { + + while (phi < 0) phi += 2*M_PI; + + //to make the math easier (and the same for each det), rotate the input phi + //BACKWARD by the phi of the det, s.t. we are looking along x-axis + phi -= center_phi; + + if (phi > M_PI) phi -= 2*M_PI; + + //then we can check easily whether it even hit the detector in phi + double det_max_phi = atan2(total_width/2,center_rho); + double det_min_phi = -det_max_phi; + + if (phi < det_min_phi || phi > det_max_phi) return std::make_pair(-1,0); + + //for theta it's not so simple, so we have to go through the typical plane-intersect method + //first thing's first: we have a fixed x for the entire detector plane: + double xhit = center_rho; + //thus we find the corresponding y and z for that fixed x, given the input theta and phi: + double yhit = xhit*tan(phi); + double zhit = sqrt(xhit*xhit+yhit*yhit)/tan(theta); + + int front_ch_hit=-1; + double ratio=0; + for (int s=0; s= front_strip_coords[s][0].GetX() && xhit <= front_strip_coords[s][0].GetX() && //Check min and max x (constant in flat) + yhit >= front_strip_coords[s][2].GetY() && yhit <= front_strip_coords[s][0].GetY() && //Check min and max y + zhit >= front_strip_coords[s][1].GetZ() && zhit <= front_strip_coords[s][0].GetZ()) //Check min and max z + { + front_ch_hit = s; + ratio = (zhit-center_z)/(length/2); + break; + } + } + + return std::make_pair(front_ch_hit,ratio); +} \ No newline at end of file diff --git a/src/ThreeStepSystem.cpp b/src/ThreeStepSystem.cpp index 33f762a..8252109 100644 --- a/src/ThreeStepSystem.cpp +++ b/src/ThreeStepSystem.cpp @@ -79,10 +79,12 @@ void ThreeStepSystem::RunSystem() { //Sample parameters double bke = generator->Gaus(m_beamDist.first, m_beamDist.second); double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second))); - double rxnPhi = 0; - double decay1Theta = GetDecayTheta(L1); + double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second); + double decay1costheta = decay1dist.GetRandomCosTheta(); + double decay1Theta = std::acos(decay1costheta); double decay1Phi = generator->Uniform(0, M_PI*2.0); - double decay2Theta = GetDecayTheta(L2); + double decay2costheta = decay2dist.GetRandomCosTheta(); + double decay2Theta = std::acos(decay2costheta); double decay2Phi = generator->Uniform(0, M_PI*2.0); double residEx = generator->Gaus(m_exDist.first, m_exDist.second); @@ -100,9 +102,22 @@ void ThreeStepSystem::RunSystem() { step1.Calculate(); step2.SetTarget(step1.GetResidual()); + if(decay1costheta == -10) { + step2.ResetEjectile(); + step2.ResetResidual(); + step3.ResetTarget(); + step3.ResetEjectile(); + step3.ResetResidual(); + return; + } step2.Calculate(); step3.SetTarget(step2.GetResidual()); + if(decay2costheta == -10) { + step3.ResetEjectile(); + step3.ResetResidual(); + return; + } step3.TurnOnResidualEloss(); step3.Calculate(); diff --git a/src/TwoStepSystem.cpp b/src/TwoStepSystem.cpp index eae927e..e7c17d2 100644 --- a/src/TwoStepSystem.cpp +++ b/src/TwoStepSystem.cpp @@ -70,8 +70,9 @@ void TwoStepSystem::RunSystem() { //Sample parameters double bke = generator->Gaus(m_beamDist.first, m_beamDist.second); double rxnTheta = acos(generator->Uniform(cos(m_theta1Range.first), cos(m_theta1Range.second))); - double rxnPhi = 0; - double decay1Theta = GetDecayTheta(L1); + double rxnPhi = generator->Uniform(m_phi1Range.first, m_phi1Range.second); + double decay1costheta = decay1dist.GetRandomCosTheta(); + double decay1Theta = std::acos(decay1costheta); double decay1Phi = generator->Uniform(0, M_PI*2.0); double residEx = generator->Gaus(m_exDist.first, m_exDist.second); @@ -86,6 +87,11 @@ void TwoStepSystem::RunSystem() { step1.Calculate(); step2.SetTarget(step1.GetResidual()); + if(decay1costheta == -10) { + step2.ResetEjectile(); + step2.ResetResidual(); + return; + } step2.TurnOnResidualEloss(); step2.Calculate(); diff --git a/src/main.cpp b/src/main.cpp index b519817..197e8af 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -2,15 +2,19 @@ #include #include "Kinematics.h" #include "SabreEfficiency.h" -#include "SabreDetector.h" +#include "AnasenEfficiency.h" #include "KinematicsExceptions.h" +#include +#include + int main(int argc, char** argv) { if(argc<2) { std::cerr<<"Incorrect number of arguments!"< detectors; - const double INNER_R = 0.0326; - const double OUTER_R = 0.1351; - const double TILT = 40.0; - const double DIST_2_TARG = -0.1245; - const double PHI_COVERAGE = 54.4; //delta phi for each det - const double PHI0 = 234.0; //center phi values for each det in array - const double PHI1 = 162.0; //# is equal to detID in channel map - const double PHI2 = 306.0; - const double PHI3 = 18.0; - const double PHI4 = 90.0; - const double DEG2RAD = M_PI/180.0; - detectors.reserve(5); - detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI0*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); - detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI1*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); - detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI2*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); - detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI3*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); - detectors.emplace_back(INNER_R,OUTER_R,PHI_COVERAGE*DEG2RAD,PHI4*DEG2RAD,TILT*DEG2RAD,DIST_2_TARG); + AnasenEfficiency anasen; + anasen.SetReactionType(calculator.GetReactionType()); + anasen.CalculateEfficiency(calculator.GetOutputName()); + //std::cout<<"Running consistency check(1=success): "<