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): "<