diff --git a/README.md b/README.md index 5600b7a..d7879dd 100644 --- a/README.md +++ b/README.md @@ -19,11 +19,7 @@ By default MASK is capable of simulating reactions of up to three steps. Here is 2. A reaction of type 2 is a reaction followed by a subsequent decay of the residual nucleus. Again, all sampling is allowed. 3. A reaction of type 3 is a reaction followed by a subsequent decay of the residual, followed by a decay of one of the products. Again, all sampling is allowed -Note that for type 2 and type 3, the decays are assumed isotropic in the center-of-mass frame of the decay. The input file requires that the user include target information, -which will be used to calculate energy loss for all of the reactants and reaction products. The target can contain layers, and each layer can be composed of a compound of elements -with a given stoichiometry. If the user wishes to not include energy loss in the kinematics, simply give all target layers a thickness of 0. Note that more layers and more thickness = more -time spent calculating energy loss. These energy loss methods are only applicable for solid targets, and should not be applied to gas or liquid targets. Energy loss calculations have a -stated uncertainty of approximately five percent. +For decays, a specific angular momentum L can be assumed. These are given in the input file as Decay1_AngularMomentum, and Decay2_AngularMomentum. This essentially modifies the center-of-mass angular distribution (as well as the lab frame). It is assumed that the decays in the center-of-mass frame are isotropic in phi (i.e. m=0). Decay1 corresponds to the first decay, if there are multiple steps, Decay2 to the second. If there are no decays, these parameters are not used (or if only one decay, Decay2_AngularMomentum is not used). The input file requires that the user include target information, which will be used to calculate energy loss for all of the reactants and reaction products. The target can contain layers, and each layer can be composed of a compound of elements with a given stoichiometry. If the user wishes to not include energy loss in the kinematics, simply give all target layers a thickness of 0. Note that more layers and more thickness = more time spent calculating energy loss. These energy loss methods are only applicable for solid targets, and should not be applied to gas or liquid targets. Energy loss calculations have a stated uncertainty of approximately five percent. The default MASK program includes a calculation of SABRE efficiency, whose methods are contained in the SabreEfficiency and SabreDetector classes. This can be disabled by modifying the main.cpp file appropriately. @@ -37,5 +33,4 @@ To run MASK simply do the following from the MASK directory: Input.txt can be replaced by any text file with the correct format. ## Requirements -MASK requires that ROOT is installed for data writting and visualization, as well as for random number generation. Testing has been done only on ROOT 6. Mileage on all other ROOT versions -will vary. \ No newline at end of file +MASK requires that ROOT is installed for data writting and visualization, as well as for random number generation. It also requires gsl to calculate Legendre Polynomials. Testing has been done only on ROOT 6. Mileage on all other ROOT versions will vary. \ No newline at end of file diff --git a/etc/DeadChannels.txt b/etc/DeadChannels.txt new file mode 100644 index 0000000..516e1ca --- /dev/null +++ b/etc/DeadChannels.txt @@ -0,0 +1,4 @@ +DetectorID RING/WEDGE Channel +4 WEDGE 4 +3 RING 2 +0 RING 13 \ No newline at end of file diff --git a/include/DeadChannelMap.h b/include/DeadChannelMap.h new file mode 100644 index 0000000..2a36426 --- /dev/null +++ b/include/DeadChannelMap.h @@ -0,0 +1,33 @@ +#ifndef DEADCHANNELMAP_H +#define DEADCHANNELMAP_H + +#include +#include +#include + +class DeadChannelMap { +public: + DeadChannelMap(); + DeadChannelMap(std::string& name); + ~DeadChannelMap(); + void LoadMapfile(std::string& name); + bool IsDead(int det, int channel, int ringwedgeFlag); + bool IsValid() { return validFlag; }; + +private: + std::unordered_map dcMap; + bool validFlag; + static constexpr int RING = 0; + static constexpr int WEDGE = 1; + + /* + Channel identifier calculated like detector*24 + 16*(RING or WEDGE) + channel + Example: + Detector 1 ring 15 + 1*24 + 16*0 + 15 = 39 + Detector 1 wedge 0 + 1*24 + 16*1 +0= 40 + */ +}; + +#endif \ No newline at end of file diff --git a/include/Kinematics.h b/include/Kinematics.h index bb26f11..24dc017 100644 --- a/include/Kinematics.h +++ b/include/Kinematics.h @@ -17,6 +17,7 @@ struct NucData { double Ex = -1; double p = -1; double theta = -1; + double theta_cm = -1; double phi = -1; int Z = -1; int A = -1; diff --git a/include/Nucleus.h b/include/Nucleus.h index 2fd0c64..8a42f52 100644 --- a/include/Nucleus.h +++ b/include/Nucleus.h @@ -21,15 +21,18 @@ public: Nucleus(int Z, int A, double px, double py, double pz, double E); ~Nucleus(); bool SetIsotope(int Z, int A); + inline void SetThetaCM(double tcm) { m_theta_cm = tcm; }; //save theta in rxn CM frame inline int GetZ() const { return m_z; }; inline int GetA() const { return m_a; }; inline double GetExcitationEnergy() const { return GetInvMass() - m_gs_mass; }; inline double GetGroundStateMass() const { return m_gs_mass; }; inline const char* GetIsotopicSymbol() const { return m_symbol.c_str(); }; + inline double GetThetaCM() const { return m_theta_cm; }; inline Nucleus& operator=(const Nucleus& rhs) { SetIsotope(rhs.GetZ(), rhs.GetA()); SetVectorCartesian(rhs.GetPx(), rhs.GetPy(), rhs.GetPz(), rhs.GetE()); + SetThetaCM(rhs.GetThetaCM()); return *this; }; @@ -46,6 +49,7 @@ public: private: int m_z, m_a; double m_gs_mass; + double m_theta_cm; std::string m_symbol; }; diff --git a/include/ReactionSystem.h b/include/ReactionSystem.h index 63fef79..dde2105 100644 --- a/include/ReactionSystem.h +++ b/include/ReactionSystem.h @@ -29,6 +29,11 @@ public: inline void SetBeamDistro(double mean, double sigma) { m_beamDist = std::make_pair(mean, sigma); }; inline void SetTheta1Range(double min, double max) { m_theta1Range = std::make_pair(min*deg2rad, max*deg2rad); }; inline void SetExcitationDistro(double mean, double sigma) { m_exDist = std::make_pair(mean, sigma); }; + inline void SetDecay1AngularMomentum(double l) { L1 = l; }; + inline void SetDecay2AngularMomentum(double l) { L2 = l; }; + + /*Sampling over angular distribution*/ + double GetDecayTheta(int L); virtual void RunSystem(); @@ -51,6 +56,7 @@ protected: 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 31a7cb2..2a3b4ae 100644 --- a/include/SabreEfficiency.h +++ b/include/SabreEfficiency.h @@ -3,12 +3,14 @@ #include "SabreDetector.h" #include "Target.h" +#include "DeadChannelMap.h" class SabreEfficiency { public: SabreEfficiency(); ~SabreEfficiency(); inline void SetReactionType(int t) { m_rxn_type = t; }; + void SetDeadChannelMap(std::string& filename) { dmap.LoadMapfile(filename); }; void CalculateEfficiency(const char* file); private: @@ -21,13 +23,13 @@ private: std::vector ringxs, ringys, ringzs; std::vector wedgexs, wedgeys, wedgezs; Target deadlayer; + DeadChannelMap dmap; //Sabre constants const double INNER_R = 0.0326; const double OUTER_R = 0.1351; const double TILT = 40.0; - //const double DIST_2_TARG = 0.14549; 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 diff --git a/input.txt b/input.txt index 54a4414..881c115 100644 --- a/input.txt +++ b/input.txt @@ -1,32 +1,28 @@ ----------Data Information---------- -OutputFile: /media/gordon/b6414c35-ec1f-4fc1-83bc-a6b68ca4325a/gwm17/test_newkine.root +OutputFile: /media/gordon/b6414c35-ec1f-4fc1-83bc-a6b68ca4325a/gwm17/12C3Hed_13N2364keV_ps_testAngMom.root SaveTree: yes SavePlots: yes ----------Reaction Information---------- ReactionType: 2 Z A (order is target, projectile, ejectile, break1, break3) -5 10 +6 12 2 3 -2 4 1 2 +1 1 ----------Target Information---------- Name: test_targ -Layers: 2 +Layers: 1 ~Layer1 -Thickness(ug/cm^2): 10 +Thickness(ug/cm^2): 41 Z A Stoich 6 12 1 0 ~ -~Layer2 -Thickness(ug/cm^2): 80 -Z A Stoich -5 10 1 -0 -~ ----------Sampling Information---------- -NumberOfSamples: 1000000 +NumberOfSamples: 100000 BeamMeanEnergy(MeV): 24 BeamEnergySigma(MeV): 0.001 -EjectileThetaMin(deg): 20.0 EjectileThetaMax(deg): 20.0 -ResidualExMean(MeV): 16.8 ResidualExSigma(MeV): 0.038 +EjectileThetaMin(deg): 3.0 EjectileThetaMax(deg): 3.0 +ResidualExMean(MeV): 3.546 ResidualExSigma(MeV): 0.02 +Decay1_AngularMomentum: 2 +Decay2_AngularMomentum: 0 -------------------------------------- diff --git a/makefile b/makefile index 1fbc246..11edf5c 100644 --- a/makefile +++ b/makefile @@ -1,8 +1,8 @@ CC=g++ ROOTGEN=rootcint -CFLAGS=-std=c++11 -g -Wall `root-config --cflags` +CFLAGS=-g -Wall `root-config --cflags` CPPFLAGS=-I ./include -LDFLAGS=`root-config --glibs` +LDFLAGS=`root-config --glibs` -lgsl ROOTINCLDIR=./ INCLDIR=./include diff --git a/src/DeadChannelMap.cpp b/src/DeadChannelMap.cpp new file mode 100644 index 0000000..8ba287e --- /dev/null +++ b/src/DeadChannelMap.cpp @@ -0,0 +1,68 @@ +#include "DeadChannelMap.h" +#include + +DeadChannelMap::DeadChannelMap() : + validFlag(false) +{ + int maxchan = 5*24+1*16+7; + for(int i=0; i>detID) { + input>>rw>>channel; + if(rw == "RING") { + this_channel = detID*24+RING*16+channel; + dcMap[this_channel] = true; + } else if(rw == "WEDGE") { + this_channel = detID*24+WEDGE*16+channel; + dcMap[this_channel] = true; + } else { + std::cerr<<"Invalid ring/wedge designation at DeadChannelMap"<second; + } else { + std::cerr<<"Error at DeadChannelMap IsDead(), invalid channel: "<>junk; } - double par1, par2; + double par1, par2, L1, L2; getline(input, junk); getline(input, junk); @@ -120,9 +120,14 @@ bool Kinematics::LoadConfig(const char* filename) { sys->SetTheta1Range(par1, par2); input>>junk>>par1>>junk>>par2; sys->SetExcitationDistro(par1, par2); + input>>junk>>L1; + input>>junk>>L2; + sys->SetDecay1AngularMomentum(L1); + sys->SetDecay2AngularMomentum(L2); sys->SetRandomGenerator(global_generator); std::cout<<"Reaction equation: "< namespace Mask { @@ -54,6 +55,22 @@ 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(gsl_sf_legendre_Pl(L, costheta), 2.0); + } while(check > probability); + + return std::acos(costheta); +} + void ReactionSystem::RunSystem() { if(!gen_set_flag) return; @@ -63,7 +80,7 @@ void ReactionSystem::RunSystem() { } if(step1.IsDecay()) { - double rxnTheta = acos(generator->Uniform(-1, 1)); + double rxnTheta = GetDecayTheta(L1); double rxnPhi = generator->Uniform(0, 2.0*M_PI); step1.SetPolarRxnAngle(rxnTheta); step1.SetAzimRxnAngle(rxnPhi); diff --git a/src/SabreEfficiency.cpp b/src/SabreEfficiency.cpp index fe42471..adfb1cf 100644 --- a/src/SabreEfficiency.cpp +++ b/src/SabreEfficiency.cpp @@ -8,6 +8,8 @@ #include #include #include +#include +#include SabreEfficiency::SabreEfficiency() : m_rxn_type(-1), deadlayer(DEADLAYER_THIN) @@ -54,6 +56,13 @@ void SabreEfficiency::CalculateEfficiency(const char* file) { std::cout<<"Loading in output from kinematics simulation: "<GetEntries(); i++) { if(++count == percent5) {//Show progress every 5% @@ -109,24 +119,38 @@ void SabreEfficiency::RunDecay(const char* file) { tree->GetEntry(i); if(eject->KE >= ENERGY_THRESHOLD) { - for(auto& det : detectors) { + 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) { - coordinates = det.GetTrajectoryCoordinates(eject->theta, eject->phi); - eject_xs.push_back(coordinates.GetX()); - eject_ys.push_back(coordinates.GetY()); - eject_zs.push_back(coordinates.GetZ()); + 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(auto& det : detectors) { - if(det.GetTrajectoryCoordinates(resid->theta, resid->phi).GetX() != 0) { - resid_xs.push_back(coordinates.GetX()); - resid_ys.push_back(coordinates.GetY()); - resid_zs.push_back(coordinates.GetZ()); + 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; } } @@ -189,15 +213,16 @@ void SabreEfficiency::Run2Step(const char* file) { std::vector b1_phis, b2_phis; std::vector b1_kes, b2_kes; - double avg_ke_per_pixel[640] = {0}; - int hits_per_pixel[640] = {0}; - //Progress tracking int percent5 = nevents*0.05; int count = 0; int npercent = 0; + int cnt_03=0, cnt_47=0, cnt_811=0, cnt_1215=0; + double avg_theta03=0, avg_theta47=0, avg_theta811=0, avg_theta1215=0; + Mask::Vec3 coords; + double thetaIncident, eloss; for(int i=0; iGetEntries(); i++) { if(++count == percent5) {//Show progress every 5% npercent++; @@ -212,14 +237,27 @@ void SabreEfficiency::Run2Step(const char* file) { 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; + + if(chan.first >= 0 && chan.first <4) { + cnt_03++; + avg_theta03 += break1->theta_cm; + } else if(chan.first >= 4 && chan.first <8) { + cnt_47++; + avg_theta47 += break1->theta_cm; + } else if(chan.first >= 8 && chan.first <12) { + cnt_811++; + avg_theta811 += break1->theta_cm; + } else if(chan.first >= 12 && chan.first <16) { + cnt_1215++; + avg_theta1215 += break1->theta_cm; + } + coords = det.GetTrajectoryCoordinates(break1->theta, break1->phi); - double thetaIncident = std::acos(coords.Dot(det.GetNormTilted())/(coords.GetR())); - double eloss = deadlayer.getEnergyLossTotal(break1->Z, break1->A, break1->KE, M_PI - thetaIncident); + 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; - int pixel = (chan.first + 16*chan.second) + 128*j; //calc pixel - avg_ke_per_pixel[pixel] += (break1->KE - eloss); - hits_per_pixel[pixel]++; b1_thetas.push_back(break1->theta); b1_phis.push_back(break1->phi); @@ -230,8 +268,16 @@ void SabreEfficiency::Run2Step(const char* file) { } if(break2->KE > ENERGY_THRESHOLD) { - for(auto& det : detectors) { - if(det.GetTrajectoryCoordinates(break2->theta, break2->phi).GetX() != 0) { + 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); @@ -239,24 +285,39 @@ void SabreEfficiency::Run2Step(const char* file) { } } } - } double b1eff = ((double) b1_thetas.size())/nevents; double b2eff = ((double) b2_thetas.size())/nevents; + double b1eff_ring03 = cnt_03/nevents; + double b1eff_ring47 = cnt_47/nevents; + double b1eff_ring811 = cnt_811/nevents; + double b1eff_ring1215 = cnt_1215/nevents; + avg_theta03 /= cnt_03; + avg_theta47 /= cnt_47; + avg_theta811 /= cnt_811; + avg_theta1215 /= cnt_1215; TParameter break1_eff("Light Breakup Efficiency", b1eff); + TParameter break1_eff_ring03("Light Breakup Efficiency Rings03", b1eff_ring03); + TParameter break1_eff_ring47("Light Breakup Efficiency Rings47", b1eff_ring47); + TParameter break1_eff_ring811("Light Breakup Efficiency Rings811", b1eff_ring811); + TParameter break1_eff_ring1215("Light Breakup Efficiency Rings1215", b1eff_ring1215); + TParameter avg_theta_cm_ring03("Avg. ThetaCM Rings03", avg_theta03); + TParameter avg_theta_cm_ring47("Avg. ThetaCM Rings47", avg_theta47); + TParameter avg_theta_cm_ring811("Avg. ThetaCM Rings811", avg_theta811); + TParameter avg_theta_cm_ring1215("Avg. ThetaCM Rings1215", avg_theta1215); TParameter break2_eff("Heavy Breakup Efficiency", b2eff); - std::ofstream output("/data1/gwm17/test_dead_pixels.txt"); - output<<"Average particle kinetic energy (MeV) per pixel (pixel = (ringch + wedgech*16) + 128*detID)"<cd(); break1_eff.Write(); + break1_eff_ring03.Write(); + break1_eff_ring47.Write(); + break1_eff_ring811.Write(); + break1_eff_ring1215.Write(); + avg_theta_cm_ring03.Write(); + avg_theta_cm_ring47.Write(); + avg_theta_cm_ring811.Write(); + avg_theta_cm_ring1215.Write(); break2_eff.Write(); input->Close(); @@ -285,6 +346,15 @@ void SabreEfficiency::Run3Step(const char* file) { int count = 0; int npercent = 0; + bool break1_flag, break3_flag, break4_flag; + 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% npercent++; @@ -292,14 +362,26 @@ void SabreEfficiency::Run3Step(const char* file) { std::cout<<"\rPercent completed: "<GetEntry(i); if(break1->KE > ENERGY_THRESHOLD) { - for(auto& det : detectors) { - if(det.GetTrajectoryCoordinates(break1->theta, break1->phi).GetX() != 0) { + 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; } } @@ -307,38 +389,88 @@ void SabreEfficiency::Run3Step(const char* file) { if(break3->KE > ENERGY_THRESHOLD) { - for(auto& det : detectors) { - if(det.GetTrajectoryCoordinates(break3->theta, break3->phi).GetX() != 0) { + 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(break4->KE > ENERGY_THRESHOLD) { - for(auto& det : detectors) { - if(det.GetTrajectoryCoordinates(break4->theta, break4->phi).GetX() != 0) { + 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(break1_flag && break3_flag && break4_flag) { + 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) { + b1b3_count++; + } else if(break1_flag && break4_flag) { + b1b4_count++; + } else if(break3_flag && break4_flag) { + 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); + } } double b1eff = ((double) b1_thetas.size())/nevents; double b3eff = ((double) b3_thetas.size())/nevents; double b4eff = ((double) b4_thetas.size())/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 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(); + b3b4_thetas->Write(); + b3b4_phis->Write(); + b3b4_kes->Write(); input->Close(); } diff --git a/src/ThreeStepSystem.cpp b/src/ThreeStepSystem.cpp index 4285ab9..33f762a 100644 --- a/src/ThreeStepSystem.cpp +++ b/src/ThreeStepSystem.cpp @@ -80,9 +80,9 @@ void ThreeStepSystem::RunSystem() { 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 = acos(generator->Uniform(-1, 1)); + double decay1Theta = GetDecayTheta(L1); double decay1Phi = generator->Uniform(0, M_PI*2.0); - double decay2Theta = acos(generator->Uniform(-1,1)); + double decay2Theta = GetDecayTheta(L2); double decay2Phi = generator->Uniform(0, M_PI*2.0); double residEx = generator->Gaus(m_exDist.first, m_exDist.second); diff --git a/src/TwoStepSystem.cpp b/src/TwoStepSystem.cpp index 54fbd43..eae927e 100644 --- a/src/TwoStepSystem.cpp +++ b/src/TwoStepSystem.cpp @@ -71,7 +71,7 @@ void TwoStepSystem::RunSystem() { 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 = acos(generator->Uniform(-1, 1)); + double decay1Theta = GetDecayTheta(L1); double decay1Phi = generator->Uniform(0, M_PI*2.0); double residEx = generator->Gaus(m_exDist.first, m_exDist.second); diff --git a/src/main.cpp b/src/main.cpp index cf461a2..b519817 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,4 +1,5 @@ #include +#include #include "Kinematics.h" #include "SabreEfficiency.h" #include "SabreDetector.h" @@ -22,7 +23,9 @@ int main(int argc, char** argv) { return 1; } SabreEfficiency sabre; + std::string mapfile = "./etc/DeadChannels.txt"; sabre.SetReactionType(calculator.GetReactionType()); + sabre.SetDeadChannelMap(mapfile); sabre.CalculateEfficiency(calculator.GetOutputName()); /*std::vector detectors;