diff --git a/Inputs/EventGenerator/sofia_238U.reaction b/Inputs/EventGenerator/sofia_238U.reaction
index 91c30af903fde508f1955d75490d4f4e10e0eedf..ff77058ce8c4af30d5e53787586a97261ae266ce 100644
--- a/Inputs/EventGenerator/sofia_238U.reaction
+++ b/Inputs/EventGenerator/sofia_238U.reaction
@@ -5,9 +5,9 @@
Beam
Particle= 238U
Energy= 135660
- ExcitationEnergy= 15 MeV
+ ExcitationEnergy= 0 MeV
SigmaEnergy= 2
- SigmaX= 6 mm
+ SigmaX= 5 mm
SigmaY= 5 mm
SigmaThetaX= 0
SigmaPhiY= 0
@@ -16,7 +16,7 @@ Beam
MeanX= 0.0 mm
MeanY= 0.0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%TwoBodyReaction
+TwoBodyReaction
Beam= 238U
Target= 208Pb
Light= 208Pb
diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
index c9ed1d005f7230999ea445726819367dc2ca509d..a937d2314309b967da2cb6c1ac170f147e4987fa 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.cxx
@@ -52,19 +52,24 @@ void TSofFissionFragment::Clear() {
fFF_TOF.clear();
fFF_Gamma.clear();
fFF_Brho.clear();
+ fFF_Rho.clear();
+ fFF_Omega.clear();
fFF_DT.clear();
fFF_Section.clear();
fFF_ThetaIn.clear();
fFF_ThetaOut.clear();
fFF_TofPosX.clear();
fFF_TofPosY.clear();
- fFF_PosX1.clear();
- fFF_PosX2.clear();
- fFF_PosX3.clear();
- fFF_PosY1.clear();
- fFF_PosY2.clear();
- fFF_PosY3.clear();
fFF_FlightPath.clear();
+ fFF_Leff.clear();
+ fFF_XB.clear();
+ fFF_XC.clear();
+ fFF_XD.clear();
+ fFF_ZB.clear();
+ fFF_ZC.clear();
+ fFF_ZD.clear();
+ fFF_X3lab.clear();
+ fFF_Z3lab.clear();
fFF_Zsum = -1;
fFF_iZsum = -1;
diff --git a/NPLib/Detectors/Sofia/TSofFissionFragment.h b/NPLib/Detectors/Sofia/TSofFissionFragment.h
index 2290d8472bb95c3fec935e723fd704de849a856e..0b30ff4a19ca7bbc91ac577a16a1fb3aa611a4b4 100644
--- a/NPLib/Detectors/Sofia/TSofFissionFragment.h
+++ b/NPLib/Detectors/Sofia/TSofFissionFragment.h
@@ -43,19 +43,24 @@ class TSofFissionFragment : public TObject {
vector<double> fFF_TOF;
vector<double> fFF_Gamma;
vector<double> fFF_Brho;
+ vector<double> fFF_Rho;
+ vector<double> fFF_Omega;
vector<double> fFF_DT;
vector<int> fFF_Section;
vector<double> fFF_ThetaIn;
vector<double> fFF_ThetaOut;
vector<double> fFF_TofPosX;
vector<double> fFF_TofPosY;
- vector<double> fFF_PosX1;
- vector<double> fFF_PosX2;
- vector<double> fFF_PosX3;
- vector<double> fFF_PosY1;
- vector<double> fFF_PosY2;
- vector<double> fFF_PosY3;
+ vector<double> fFF_XB;
+ vector<double> fFF_XC;
+ vector<double> fFF_XD;
+ vector<double> fFF_ZB;
+ vector<double> fFF_ZC;
+ vector<double> fFF_ZD;
+ vector<double> fFF_X3lab;
+ vector<double> fFF_Z3lab;
vector<double> fFF_FlightPath;
+ vector<double> fFF_Leff;
double fFF_Zsum;
int fFF_iZsum;
@@ -92,27 +97,32 @@ class TSofFissionFragment : public TObject {
inline void SetTOF(double val){fFF_TOF.push_back(val);};//!
inline void SetGamma(double val){fFF_Gamma.push_back(val);};//!
inline void SetBrho(double val){fFF_Brho.push_back(val);};//!
+ inline void SetRho(double val){fFF_Rho.push_back(val);};//!
+ inline void SetOmega(double val){fFF_Omega.push_back(val);};//!
inline void SetDT(double val){fFF_DT.push_back(val);};//!
inline void SetSection(int val){fFF_Section.push_back(val);};//!
inline void SetThetaIn(double val){fFF_ThetaIn.push_back(val);};//!
inline void SetThetaOut(double val){fFF_ThetaOut.push_back(val);};//!
inline void SetTofPosX(double val){fFF_TofPosX.push_back(val);};//!
inline void SetTofPosY(double val){fFF_TofPosY.push_back(val);};//!
- inline void SetPosX1(double val){fFF_PosX1.push_back(val);};//!
- inline void SetPosX2(double val){fFF_PosX2.push_back(val);};//!
- inline void SetPosX3(double val){fFF_PosX3.push_back(val);};//!
- inline void SetPosY1(double val){fFF_PosY1.push_back(val);};//!
- inline void SetPosY2(double val){fFF_PosY2.push_back(val);};//!
- inline void SetPosY3(double val){fFF_PosY3.push_back(val);};//!
inline void SetFlightPath(double val){fFF_FlightPath.push_back(val);};//!
+ inline void SetLeff(double val){fFF_Leff.push_back(val);};//!
+ inline void SetPosXB(double val){fFF_XB.push_back(val);};//!
+ inline void SetPosXC(double val){fFF_XC.push_back(val);};//!
+ inline void SetPosXD(double val){fFF_XD.push_back(val);};//!
+ inline void SetPosZB(double val){fFF_ZB.push_back(val);};//!
+ inline void SetPosZC(double val){fFF_ZC.push_back(val);};//!
+ inline void SetPosZD(double val){fFF_ZD.push_back(val);};//!
+ inline void SetPosX3lab(double val){fFF_X3lab.push_back(val);};//!
+ inline void SetPosZ3lab(double val){fFF_Z3lab.push_back(val);};//!
////////////////////// GETTERS ////////////////////////
int GetMult() {return fFF_Z.size();}//!
int GetMultTofPos() {return fFF_TofPosY.size();}//!
- int GetMultMwpc1() {return fFF_PosY1.size();}//!
+ /*int GetMultMwpc1() {return fFF_PosY1.size();}//!
int GetMultMwpc2() {return fFF_PosY2.size();}//!
- int GetMultMwpc3() {return fFF_PosY3.size();}//!
+ int GetMultMwpc3() {return fFF_PosY3.size();}//!*/
inline double GetZsum() const {return fFF_Zsum;}//!
inline int GetiZsum() const {return fFF_iZsum;}//!
inline int GetPlastic(int i) const {return fFF_Plastic[i];}//!
@@ -124,19 +134,24 @@ class TSofFissionFragment : public TObject {
inline double GetTOF(int i) const {return fFF_TOF[i];}//!
inline double GetGamma(int i) const {return fFF_Gamma[i];}//!
inline double GetBrho(int i) const {return fFF_Brho[i];}//!
+ inline double GetRho(int i) const {return fFF_Rho[i];}//!
+ inline double GetOmega(int i) const {return fFF_Omega[i];}//!
inline double GetDT(int i) const {return fFF_DT[i];}//!
inline double GetSection(int i) const {return fFF_Section[i];}//!
inline double GetThetaIn(int i) const {return fFF_ThetaIn[i];}//!
inline double GetThetaOut(int i) const {return fFF_ThetaOut[i];}//!
inline double GetTofPosX(int i) const {return fFF_TofPosX[i];}//!
inline double GetTofPosY(int i) const {return fFF_TofPosY[i];}//!
- inline double GetPosX1(int i) const {return fFF_PosX1[i];}//!
- inline double GetPosX2(int i) const {return fFF_PosX2[i];}//!
- inline double GetPosX3(int i) const {return fFF_PosX3[i];}//!
- inline double GetPosY1(int i) const {return fFF_PosY1[i];}//!
- inline double GetPosY2(int i) const {return fFF_PosY2[i];}//!
- inline double GetPosY3(int i) const {return fFF_PosY3[i];}//!
inline double GetFlightPath(int i) const {return fFF_FlightPath[i];}//!
+ inline double GetLeff(int i) const {return fFF_Leff[i];}//!
+ inline double GetPosXB(int i) const {return fFF_XB[i];}//!
+ inline double GetPosXC(int i) const {return fFF_XC[i];}//!
+ inline double GetPosXD(int i) const {return fFF_XD[i];}//!
+ inline double GetPosZB(int i) const {return fFF_ZB[i];}//!
+ inline double GetPosZC(int i) const {return fFF_ZC[i];}//!
+ inline double GetPosZD(int i) const {return fFF_ZD[i];}//!
+ inline double GetPosX3lab(int i) const {return fFF_X3lab[i];}//!
+ inline double GetPosZ3lab(int i) const {return fFF_Z3lab[i];}//!
//////////////////////////////////////////////////////////////
diff --git a/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx b/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx
index f5f1a5e66097486ff7b81b4a3415477cfcc7411f..a1a69269ef1f935bb44a17342a45c6454da14a58 100644
--- a/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx
+++ b/NPLib/Detectors/Sofia/TSofTwimPhysics.cxx
@@ -193,7 +193,7 @@ void TSofTwimPhysics::BuildPhysicalEvent() {
double p1 = (12*xysum -xsum*ysum)/(12*x2sum - xsum*xsum);
Theta1 = atan(p1);
- Theta.push_back(-Theta1);
+ Theta.push_back(Theta1);
}
if(Esec2>0){
Esec2 = Cal->ApplyCalibration("SofTwim/SEC2_ALIGN",Esec2);
@@ -217,7 +217,7 @@ void TSofTwimPhysics::BuildPhysicalEvent() {
double p1 = (12*xysum -xsum*ysum)/(12*x2sum - xsum*xsum);
Theta2 = atan(p1);
- Theta.push_back(-Theta2);
+ Theta.push_back(Theta2);
}
if(Esec3>0){
Esec3 = Cal->ApplyCalibration("SofTwim/SEC3_ALIGN",Esec3);
@@ -241,7 +241,7 @@ void TSofTwimPhysics::BuildPhysicalEvent() {
double p1 = (12*xysum -xsum*ysum)/(12*x2sum - xsum*xsum);
Theta3 = atan(p1);
- Theta.push_back(-Theta3);
+ Theta.push_back(Theta3);
}
if(Esec4>0){
Esec4 = Cal->ApplyCalibration("SofTwim/SEC4_ALIGN",Esec4);
@@ -265,7 +265,7 @@ void TSofTwimPhysics::BuildPhysicalEvent() {
double p1 = (12*xysum -xsum*ysum)/(12*x2sum - xsum*xsum);
Theta4 = atan(p1);
- Theta.push_back(-Theta4);
+ Theta.push_back(Theta4);
}
m_Beta = -1;
diff --git a/NPSimulation/Process/FissionDecay.cc b/NPSimulation/Process/FissionDecay.cc
index e126a46cb95975954baa6d3bd92e00cc688cadeb..e8776516478b87672f84c0ddf31f7f3b0e8a6fab 100644
--- a/NPSimulation/Process/FissionDecay.cc
+++ b/NPSimulation/Process/FissionDecay.cc
@@ -183,12 +183,12 @@ void FissionDecay::DoIt(const G4FastTrack& fastTrack,G4FastStep& fastStep){
if(Zsum == m_CompoundParticle.GetZ())
m_FissionConditions->SetNeutronMultiplicity(m_CompoundParticle.GetA()-Asum);
- for(unsigned int i = 0 ; i < size ; i++){
+ for(int i=size-1; i>=0; --i){
// Get the decaying particle
int FFZ = FissionFragment[i].GetZ();
int FFA = FissionFragment[i].GetA();
FissionFragmentDef=NULL;
-
+
// Set the momentum direction
G4ThreeVector Momentum (DPx[i],DPy[i],DPz[i]);
Momentum=Momentum.unit();
diff --git a/Projects/Sofia/sofia.detector b/Projects/Sofia/sofia.detector
index 36c9551bc5782ed58fd03b19d0376cdb35358eb6..3cba4f8d2c222f6468809d00271428e3f84154c0 100644
--- a/Projects/Sofia/sofia.detector
+++ b/Projects/Sofia/sofia.detector
@@ -12,7 +12,7 @@ SofTofW
R= 4.5 m
THETA= -18 deg
PHI= 0 deg
- Yoffset= 30.0 mm
+ Yoffset= 0.0 mm
Build_GLAD= 1
GLAD_TiltAngle= 14 deg
GLAD_DistanceFromTarget= 3.3 m
diff --git a/Projects/e793s/Analysis.cxx b/Projects/e793s/Analysis.cxx
index c04cbfdd958e03e3fbb5f0552925b3f42c6cb1d5..aa1eabc63a3665dd96697432752483cd743fcdbb 100755
--- a/Projects/e793s/Analysis.cxx
+++ b/Projects/e793s/Analysis.cxx
@@ -227,9 +227,11 @@ void Analysis::TreatEvent(){
/**/
if(isSim && !isPhaseSpace){
- if(M2->Si_E[countMust2]>0 && M2->CsI_E[countMust2]<=0){
- /* IS a count in DSSD, but NOT in CsI to exclude punch through */
- /* May need further gating in Si_E & _T */
+ if(M2->TelescopeNumber[countMust2]<5){
+ if(M2->Si_E[countMust2]>0 && // DSSD count
+ M2->CsI_E[countMust2]<=0 && // No CsI count
+ M2->Si_T[countMust2]<460 // Triton kinematic line, not punch through
+ ){
ThetaCM_detected_MM->Fill(ReactionConditions->GetThetaCM());
ThetaLab_detected_MM->Fill(ReactionConditions->GetTheta(0));
@@ -237,6 +239,15 @@ void Analysis::TreatEvent(){
ThetaCM_detected_MMX[MMX]->Fill(ReactionConditions->GetThetaCM());
ThetaLab_detected_MMX[MMX]->Fill(ReactionConditions->GetTheta(0));
}
+ } else {
+ //No triton requirement for MM5
+ ThetaCM_detected_MM->Fill(ReactionConditions->GetThetaCM());
+ ThetaLab_detected_MM->Fill(ReactionConditions->GetTheta(0));
+
+ int MMX = TelescopeNumber-1;
+ ThetaCM_detected_MMX[MMX]->Fill(ReactionConditions->GetThetaCM());
+ ThetaLab_detected_MMX[MMX]->Fill(ReactionConditions->GetTheta(0));
+ }
}
/**/
diff --git a/Projects/e793s/Analysis.h b/Projects/e793s/Analysis.h
index 071c320aaa57eb6023ed4aa93f71ad0d8533c6e8..150d05057f51ff658658845895a4bf71bec4fe36 100755
--- a/Projects/e793s/Analysis.h
+++ b/Projects/e793s/Analysis.h
@@ -42,7 +42,7 @@
#include <TMath.h>
#include <bitset>
-int NumThetaAngleBins = 900;// 180 = 1 deg, 360 = 0.5 deg, 900 = 0.2 deg, 1800 = 0.1 deg
+int NumThetaAngleBins = 1800;// 180 = 1 deg, 360 = 0.5 deg, 900 = 0.2 deg, 1800 = 0.1 deg
auto ThetaCM_emmitted = new TH1F("ThetaCM_emmitted","ThetaCM_emmitted",NumThetaAngleBins,0,180);
auto ThetaCM_detected_MM = new TH1F("ThetaCM_detected_MM","ThetaCM_detected_MM",NumThetaAngleBins,0,180);
diff --git a/Projects/e793s/macro/CS2_dt.h b/Projects/e793s/macro/CS2_dt.h
index a479ea275bdea7a908d86d33d8c4736453ca1ba7..87614b9b695c87058152dcfd1b1090e8c67736ac 100644
--- a/Projects/e793s/macro/CS2_dt.h
+++ b/Projects/e793s/macro/CS2_dt.h
@@ -94,8 +94,9 @@ void CS_Diagnosis(){
void CS(){
/* Overload function */
cout << "- CS(stateE, stateSp, orb_l, orb_j, nodes) "<< endl;
- cout << "---- 1.945, p1/2 = CS(1.945, 1, 0, 0.5) "<< endl;
+ cout << "---- 1.945, s1/2 = CS(1.945, 1, 0, 0.5) "<< endl;
cout << "---- 1.945, d3/2 = CS(1.945, 1, 2, 1.5) "<< endl;
+ cout << "---- 0.691, f7/2 = CS(0.691, 4, 3, 3.5) "<< endl;
}
////////////////////////////////////////////////////////////////////////////////
@@ -536,23 +537,19 @@ vector<vector<double>> GetExpDiffCross(double Energy){
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
- // TEMPORARY FIX FOR EX:THETACM PROBLEM!!
- // Assuming angular bins are 2-20 in 1 degree sections...
- if(means_dt[indexE] > 4.0){
- numAngleBins-= 4; //Up to 16 deg
- }
- else if(means_dt[indexE] > 3.0){
- numAngleBins-= 5; //Up to 15 deg
- }
- else if(means_dt[indexE] > 2.0){
- numAngleBins-= 6; //Up to 14 deg
- }
- else if(means_dt[indexE] > 1.0){
- numAngleBins-= 8; //Up to 12 deg
- }
- else if(means_dt[indexE] > 0.0){
- numAngleBins-=12; //Up to 8 deg
- }
+ /* TEMPORARY FIX FOR EX:THETACM PROBLEM!! */
+ double thetaMax = 0.23*pow(means_dt[indexE],2)
+ + 0.76*means_dt[indexE]
+ + 9.01;
+ thetaMax = floor(thetaMax);
+
+ double thetaMin = 0.10*pow(means_dt[indexE],2)
+ + 0.06*means_dt[indexE]
+ + 1.58;
+ thetaMin = ceil(thetaMin);
+
+ numAngleBins = (int)thetaMax - (int)thetaMin;
+ firstAngle = (int)thetaMin;
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
diff --git a/Projects/e793s/macro/DrawPlots.h b/Projects/e793s/macro/DrawPlots.h
index 4456c7e58f5d4cd6ac88822052ea0ede20f41593..44eac574080ab87a6e67a30d3edd047ada187414 100755
--- a/Projects/e793s/macro/DrawPlots.h
+++ b/Projects/e793s/macro/DrawPlots.h
@@ -1199,6 +1199,39 @@ void ExThetaLab(){
Sn->Draw();
}
+void ExThetaCM(){
+ string gate = timegate
+ + " && " + det_gate
+ + " && Ex@.size()==1";
+ if(reactionName=="47K(d,t)"){
+ gate = gate + " && cutTritons && cutTime";
+ }
+
+ TCanvas *diagnoseTheta = new TCanvas("diagnoseTheta","diagnoseTheta",1000,1000);
+ chain->Draw(
+ "Ex:ThetaCM>>thetaHistCM(360,0,180,180,-1,8)",
+ gate.c_str(), "colz");
+ TH1F* thetaHistCM = (TH1F*) gDirectory->Get("thetaHistCM");
+ thetaHistCM->GetXaxis()->SetTitle("#theta_{CM} [deg]");
+ thetaHistCM->GetYaxis()->SetTitle("Ex [MeV]");
+ thetaHistCM->SetTitle("Theta dependance testing");
+
+ diagnoseTheta->Update();
+
+ TLine *l0000 = new TLine(100., 0.000, 160., 0.000);
+ l0000->SetLineStyle(kDashed);
+ l0000->SetLineColor(kRed);
+ l0000->Draw();
+
+ if(reactionName=="47K(d,p)"){
+ TLine *Sn = new TLine(100., 4.644, 160., 4.644);
+ Sn->SetLineStyle(kDashed);
+ Sn->SetLineColor(kRed);
+ Sn->Draw();
+ }
+}
+
+
void ExThetaLab(double gamma, double width){
TCanvas *diagnoseTheta = new TCanvas("diagnoseTheta","diagnoseTheta",1000,1000);
string gate = timegate
@@ -2864,3 +2897,120 @@ void Figure_TopGamma_BottomParticle(double gammaBinWidth, double particleBinWidt
cFig_GGSP->Update();
}
+
+void Figure_SolidAngle_dt(){
+ //string histname = "SolidAngle_CM_MM";
+ string histname = "SolidAngle_Lab_MM";
+
+ TCanvas* canv = new TCanvas("canv","canv",1000,1000);
+ gStyle->SetPadLeftMargin(0.10);
+ gStyle->SetPadRightMargin(0.03);
+ gStyle->SetOptStat(0);
+
+ TFile *file6 = new TFile("./SolidAngle_HistFiles/SAHF_18Oct22_47Kdt_6000.root","READ");
+ TFile *file5 = new TFile("./SolidAngle_HistFiles/SAHF_18Oct22_47Kdt_5000.root","READ");
+ TFile *file4 = new TFile("./SolidAngle_HistFiles/SAHF_18Oct22_47Kdt_4000.root","READ");
+ TFile *file3 = new TFile("./SolidAngle_HistFiles/SAHF_18Oct22_47Kdt_3000.root","READ");
+ TFile *file2 = new TFile("./SolidAngle_HistFiles/SAHF_18Oct22_47Kdt_2000.root","READ");
+ TFile *file1 = new TFile("./SolidAngle_HistFiles/SAHF_18Oct22_47Kdt_1000.root","READ");
+ TFile *file0 = new TFile("./SolidAngle_HistFiles/SAHF_18Oct22_47Kdt_0000.root","READ");
+
+ TH1F* h6 = (TH1F*)file6->Get(histname.c_str());
+ //h6->GetXaxis()->SetTitle("#theta_{CM}");
+ h6->GetXaxis()->SetTitle("#theta_{Lab}");
+ h6->GetYaxis()->SetTitle("Solid Angle [(sr)?]");
+ h6->SetTitle("Simulated triton solid angle variation with 46K state energy");
+ h6->GetXaxis()->SetRangeUser(0.,30.);
+ h6->GetYaxis()->SetRangeUser(0.,0.002);
+ h6->SetLineColor(kCyan+4);
+ h6->SetFillColor(kCyan+4);
+ h6->SetLineWidth(2);
+ h6->SetFillStyle(3002);
+ h6->Draw("hist");
+
+ TH1F* h5 = (TH1F*)file5->Get(histname.c_str());
+ h5->SetLineColor(kCyan-1);
+ h5->SetFillColor(kCyan-1);
+ h5->SetLineWidth(2);
+ h5->SetFillStyle(3002);
+ h5->Draw("hist same");
+
+ TH1F* h4 = (TH1F*)file4->Get(histname.c_str());
+ h4->SetLineColor(kCyan-5);
+ h4->SetFillColor(kCyan-5);
+ h4->SetLineWidth(2);
+ h4->SetFillStyle(3002);
+ h4->Draw("hist same");
+
+ TH1F* h3 = (TH1F*)file3->Get(histname.c_str());
+ h3->SetLineColor(kCyan-8);
+ h3->SetFillColor(kCyan-8);
+ h3->SetLineWidth(2);
+ h3->SetFillStyle(3002);
+ h3->Draw("hist same");
+
+ TH1F* h2 = (TH1F*)file2->Get(histname.c_str());
+ h2->SetLineColor(kCyan-9);
+ h2->SetFillColor(kCyan-9);
+ h2->SetLineWidth(2);
+ h2->SetFillStyle(3002);
+ h2->Draw("hist same");
+
+ TH1F* h1 = (TH1F*)file1->Get(histname.c_str());
+ h1->SetLineColor(kCyan-7);
+ h1->SetFillColor(kCyan-7);
+ h1->SetLineWidth(2);
+ h1->SetFillStyle(3002);
+ h1->Draw("hist same");
+
+ TH1F* h0 = (TH1F*)file0->Get(histname.c_str());
+ h0->SetLineColor(kCyan-0);
+ h0->SetFillColor(kCyan-0);
+ h0->SetLineWidth(2);
+ h0->SetFillStyle(3002);
+ h0->Draw("hist same");
+
+ TLatex *n0 = new TLatex(.5,.5,"0.0 MeV");
+ n0->SetTextColor(kCyan-0);
+ n0->SetTextSize(0.05);
+ n0->SetX(22.);
+ n0->SetY(0.0018);
+ n0->Draw("same");
+ TLatex *n1 = new TLatex(.5,.5,"1.0 MeV");
+ n1->SetTextColor(kCyan-7);
+ n1->SetTextSize(0.05);
+ n1->SetX(22.);
+ n1->SetY(0.0017);
+ n1->Draw("same");
+ TLatex *n2 = new TLatex(.5,.5,"2.0 MeV");
+ n2->SetTextColor(kCyan-9);
+ n2->SetTextSize(0.05);
+ n2->SetX(22.);
+ n2->SetY(0.0016);
+ n2->Draw("same");
+ TLatex *n3 = new TLatex(.5,.5,"3.0 MeV");
+ n3->SetTextColor(kCyan-8);
+ n3->SetTextSize(0.05);
+ n3->SetX(22.);
+ n3->SetY(0.0015);
+ n3->Draw("same");
+ TLatex *n4 = new TLatex(.5,.5,"4.0 MeV");
+ n4->SetTextColor(kCyan-5);
+ n4->SetTextSize(0.05);
+ n4->SetX(22.);
+ n4->SetY(0.0014);
+ n4->Draw("same");
+ TLatex *n5 = new TLatex(.5,.5,"5.0 MeV");
+ n5->SetTextColor(kCyan-1);
+ n5->SetTextSize(0.05);
+ n5->SetX(22.);
+ n5->SetY(0.0013);
+ n5->Draw("same");
+ TLatex *n6 = new TLatex(.5,.5,"6.0 MeV");
+ n6->SetTextColor(kCyan+4);
+ n6->SetTextSize(0.05);
+ n6->SetX(22.);
+ n6->SetY(0.0012);
+ n6->Draw("same");
+
+}
diff --git a/Projects/e793s/macro/KnownPeakFitter.h b/Projects/e793s/macro/KnownPeakFitter.h
index 42b6a3ccd6cfe4d6ac0f5476bbf7cc417a2fba2d..f4f66e89c990396dd31411fc02525c27dca15d35 100644
--- a/Projects/e793s/macro/KnownPeakFitter.h
+++ b/Projects/e793s/macro/KnownPeakFitter.h
@@ -24,10 +24,17 @@ array<double,numPeaks> means = { 0.000,
5.82
};
-const int numPeaks_dt = 3;
+const int numPeaks_dt = 8;//9;
array<double,numPeaks_dt> means_dt = { 0.000,
+ //0.587, //from lit
+ 0.691, //from lit
+ //0.886, //from lit
+ 1.738, //from lit
1.945,
+ 2.26 , //from lit
+ 2.76 , //from lit
3.344,
+ 4.2
};
@@ -371,3 +378,9 @@ void FitKnownPeaks(TH1F* hist){
//Shell function to call Rtrn_Arry without writing vector<vector<double>> to screen
vector<vector<double>> shell = FitKnownPeaks_RtrnArry(hist,0.0);
}
+
+void FitKnownPeaks_dt(TH1F* hist){
+ //Shell function to call Rtrn_Arry without writing vector<vector<double>> to screen
+ vector<vector<double>> shell = FitKnownPeaks_dt_RtrnArry(hist,0.0);
+}
+
diff --git a/Projects/s455/Analysis.cxx b/Projects/s455/Analysis.cxx
index 421c7972339d263e74b35a7b1b11c087487bcd7e..1ada8376f1b7f451fbf612f6bd3d8ced122a2d67 100644
--- a/Projects/s455/Analysis.cxx
+++ b/Projects/s455/Analysis.cxx
@@ -37,8 +37,10 @@ struct TofPair
double tof = -1000;
double velocity = -1;
double beta = -1;
+ double gamma = -1;
double theta_in = -10;
double theta_out = -10;
+ double psi = -10;
int plastic = -1;
//
int isLorR = -1;
@@ -50,6 +52,10 @@ struct TofPair
// *** isUorD = 2 -> Down
int section = -1;
double Esec = -1;
+ double Z = 0;
+ int iZ = 0;
+ double AoQ = 0;
+ double A = 0;
double DT = -100;
double x2twim = -1000;
double x2 = -1000;
@@ -57,10 +63,18 @@ struct TofPair
double y3 = -1000;
double x3lab = 0;
double z3lab = 0;
+ double xb = 0;
double xc = 0;
+ double xd = 0;
double yc = 0;
+ double zb = 0;
double zc = 0;
+ double zd = 0;
double flight_path = 0;
+ double Leff = 0;
+ double rho = 0;
+ double Brho = 0;
+ double omega = 0;
};
@@ -540,12 +554,6 @@ void Analysis::FissionFragmentAnalysis(){
for(int i=0; i<2; i++){
int section = TofHit[i].section;
- /*double drift_time;
- if(section==1) drift_time = DT1;
- if(section==2) drift_time = DT2;
- if(section==3) drift_time = DT3;
- if(section==4) drift_time = DT4;*/
-
double DT_eval;
if(section<3) DT_eval=55;
if(section>2) DT_eval=-55;
@@ -561,164 +569,131 @@ void Analysis::FissionFragmentAnalysis(){
double Theta0 = 20.*deg;
double XA = 0;
double ZA = 2272;
- double Zfission = 335.;
- double ZGlad_entrance = 2694.;
+ //double ZGlad_entrance = 2694.;
+ double ZGlad_entrance = 2694.+540.5;
double Leff_init = 2150.;
- //double ZG = ZGlad_entrance+Leff_init/2;
- double ZG = ZGlad_entrance+540;
+ //double ZG = ZGlad_entrance+1654.;
+ double ZG = ZGlad_entrance+1113.5;
double ZMW3 = 8450;//8483;
double XMW3 = -(ZMW3-ZG)*tan(Theta0);
double ZMW2 = 2576;
double X3lab = 0;
double Z3lab = 0;;
double Tilt = 14.*deg;
- TVector3 vOut;
TVector3 vZ = TVector3(0,0,1);
- TVector3 vC;
- TVector3 InitPos[2];
- TVector3 InitDir[2];
- TVector3 FinalPos[2];
double XC = 0;
double ZC = 0;
- double Leff = 0;
+ double XB = 0;
+ double ZB = 0;
+ double XD = 0;
+ double ZD = 0;
+ double MagB = m_GladField->GetB()/1000;
for(int i=0; i<2; i++){
XA = TofHit[i].DT;
if(XA != -1e6){
+ // *** Extroplate to C position *** //
XC = (XA+(ZG-ZA)*tan(TofHit[i].theta_in)) / (1-tan(Tilt)*tan(TofHit[i].theta_in));
ZC = ZG + XC*tan(Tilt);
-
- /*int ix = (int)(XC - m_GladField->GetXmin())/50;
- int iy = (int)(ZC - m_GladField->GetYmin())/50;
- Leff = m_GladField->GetLeff(ix,iy);
- ZG = ZGlad_entrance + Leff/2;
- XC = (XA+(ZG-ZA)*tan(TofHit[i].theta_in)) / (1-tan(Tilt)*tan(TofHit[i].theta_in));
- ZC = ZG + XC*tan(Tilt);*/
TofHit[i].xc = XC;
- TofHit[i].yc = TofHit[i].y*0.5;
+ TofHit[i].yc = 0;//TofHit[i].y*0.5;
TofHit[i].zc = ZC;
+ int ix, iy;
+ ix = (int)(TofHit[0].xc - m_GladField->GetXmin())/50;
+ iy = (int)(TofHit[0].yc - m_GladField->GetYmin())/50;
+ TofHit[i].Leff = m_GladField->GetLeff(ix,iy);
+
X3lab = TofHit[i].x3*cos(Theta0) + XMW3;
Z3lab = TofHit[i].x3*sin(Theta0) + ZMW3;
TofHit[i].x3lab = X3lab;
TofHit[i].z3lab = Z3lab;
-
- InitPos[i] = TVector3(XA,0,ZA);
- InitDir[i] = TVector3(sin(TofHit[i].theta_in),0,cos(TofHit[i].theta_in));
- FinalPos[i] = TVector3(X3lab,0,Z3lab);
-
- vC = TVector3(TofHit[i].xc,0,TofHit[i].zc);
- vOut = TVector3(X3lab-TofHit[i].xc,0,Z3lab-TofHit[i].zc);
-
- double PathLength = vC.Mag() + vOut.Mag() + 74.;
+ TVector3 vC = TVector3(TofHit[i].xc,0,TofHit[i].zc);
+ TVector3 vOut = TVector3(X3lab-TofHit[i].xc,0,Z3lab-TofHit[i].zc);
+ double angle = -vZ.Angle(vOut);
+ TofHit[i].theta_out = angle;
+ TofHit[i].psi = TofHit[i].theta_in - TofHit[i].theta_out;
+ TofHit[i].rho = TofHit[i].Leff/(2*sin(0.5*TofHit[i].psi)*cos(Tilt-0.5*TofHit[i].psi));
+ TofHit[i].Brho = MagB*TofHit[i].rho;
+
+ // *** Extrapolate to B position *** //
+ double ZI = ZG - TofHit[i].Leff/(2*cos(Tilt));
+ XB = (XA+(ZI-ZA)*tan(TofHit[i].theta_in)) / (1-tan(Tilt)*tan(TofHit[i].theta_in));
+ ZB = ZI + XB*tan(Tilt);
+ TofHit[i].xb = XB;
+ TofHit[i].zb = ZB;
+ TVector3 vB = TVector3(XB,0,ZB);
+
+ // *** Extrapolate to D position *** //
+ //XD = XC + TofHit[i].Leff/2*sin(angle)/cos(angle+Tilt);
+ XD = XC + TofHit[i].Leff/2*tan(angle+Tilt)*cos(Tilt);
+ //ZD = ZC + TofHit[i].Leff/2*cos(angle)/cos(angle+Tilt);
+ ZD = ZC + TofHit[i].Leff/2*cos(Tilt);
+ TofHit[i].xd = XD;
+ TofHit[i].zd = ZD;
+ TVector3 vD = TVector3(XD,0,ZD);
+
+ TVector3 v1 = TVector3(XB,0,ZB);
+ TVector3 v3 = TVector3(X3lab-XD,0,Z3lab-ZD);
+ TofHit[i].omega = abs(2.*asin(sqrt(pow(XD-XB,2) + pow(ZD-ZB,2))/(2*TofHit[i].rho)));
+ double Path1 = v1.Mag();
+ double Path2 = TofHit[i].rho*TofHit[i].omega;
+ double Path3 = v3.Mag();
+ double PathLength = Path1 + Path2 + Path3 + 74.;
PathLength = PathLength/1000.;
- double angle = vZ.Angle(vOut);
TofHit[i].flight_path = PathLength;
TofHit[i].velocity = PathLength/TofHit[i].tof;
TofHit[i].beta = TofHit[i].velocity * m/ns / NPUNITS::c_light;
- TofHit[i].theta_out = angle;
- TofHit[i].x2twim = XA + (ZMW2-ZA)*tan(TofHit[i].theta_in);
+ TofHit[i].x2twim = XA;
+
+ double Z = TofHit[i].Esec;
+ Z = fZff_p0 + fZff_p1*Z + fZff_p2*Z*Z;
+ Z = sqrt(Z);
+ int iZ = (int) round(Z);
+ TofHit[i].Z = Z;
+ TofHit[i].iZ = iZ;
+
+ TofHit[i].gamma = 1. / sqrt(1 - TofHit[i].beta*TofHit[i].beta);
+ TofHit[i].AoQ = TofHit[i].Brho / (3.10761 * TofHit[i].beta * TofHit[i].gamma);
+ TofHit[i].A = TofHit[i].AoQ * TofHit[i].iZ;
+
}
}
- Z1 = TofHit[0].Esec;
- Z2 = TofHit[1].Esec;
-
-
- if(Z1>0 && Z2>0){
- Z1 = fZff_p0 + fZff_p1*Z1 + fZff_p2*Z1*Z1;
- Z2 = fZff_p0 + fZff_p1*Z2 + fZff_p2*Z2*Z2;
-
- Z1 = sqrt(Z1);
- Z2 = sqrt(Z2);
-
- Zsum = Z1+Z2;
+ // *** Filling the Fission Fragment Tree *** //
+ for(int i=0; i<2; i++){
+ SofFF->SetTOF(TofHit[i].tof);
+ SofFF->SetTofPosX(TofHit[i].x);
+ SofFF->SetTofPosY(TofHit[i].y);
+ SofFF->SetPlastic(TofHit[i].plastic);
+ SofFF->SetPosXB(TofHit[i].xb);
+ SofFF->SetPosXC(TofHit[i].xc);
+ SofFF->SetPosXD(TofHit[i].xd);
+ SofFF->SetPosZB(TofHit[i].zb);
+ SofFF->SetPosZC(TofHit[i].zc);
+ SofFF->SetPosZD(TofHit[i].zd);
+ SofFF->SetPosX3lab(TofHit[i].x3lab);
+ SofFF->SetPosZ3lab(TofHit[i].z3lab);
+ SofFF->SetThetaIn(TofHit[i].theta_in/deg);
+ SofFF->SetThetaOut(TofHit[i].theta_out/deg);
+ SofFF->SetBeta(TofHit[i].beta);
+ SofFF->SetGamma(TofHit[i].gamma);
+ SofFF->SetiZ(TofHit[i].iZ);
+ SofFF->SetZ(TofHit[i].Z);
+ SofFF->SetAoQ(TofHit[i].AoQ);
+ SofFF->SetA(TofHit[i].A);
+ SofFF->SetBrho(TofHit[i].Brho);
+ SofFF->SetRho(TofHit[i].rho);
+ SofFF->SetOmega(TofHit[i].omega);
+ SofFF->SetDT(TofHit[i].DT);
+ SofFF->SetSection(TofHit[i].section);
+ SofFF->SetLeff(TofHit[i].Leff);
+ SofFF->SetFlightPath(TofHit[i].flight_path);
- iZ1 = (int) round(Z1);
- iZ2 = (int) round(Z2);
- iZsum = iZ1 + iZ2;
}
-
- double MagB = m_GladField->GetB()/1000;
- int ix, iy;
- ix = (int)(TofHit[0].xc - m_GladField->GetXmin())/50;
- iy = (int)(TofHit[0].yc - m_GladField->GetYmin())/50;
- double Leff1 = m_GladField->GetLeff(ix,iy);
- ix = (int)(TofHit[1].xc - m_GladField->GetXmin())/50;
- iy = (int)(TofHit[1].yc - m_GladField->GetYmin())/50;
- double Leff2 = m_GladField->GetLeff(ix,iy);
- double rho1 = Leff1/abs(2*sin(0.5*(TofHit[0].theta_out-TofHit[0].theta_in))*cos(Tilt-0.5*(TofHit[0].theta_out-TofHit[0].theta_in)));
- double rho2 = Leff2/abs(2*sin(0.5*(TofHit[1].theta_out-TofHit[1].theta_in))*cos(Tilt-0.5*(TofHit[1].theta_out-TofHit[1].theta_in)));
- double Brho1 = MagB*rho1;
- double Brho2 = MagB*rho2;
-
- /*double Brho1 = 0;
- double Brho2 = 0;
- Brho1 = m_GladField->FindBrho(InitPos[0], InitDir[0], FinalPos[0]);
- Brho2 = m_GladField->FindBrho(InitPos[1], InitDir[1], FinalPos[1]);
- //cout << Brho1 << " " << Brho2 << endl;*/
-
- Beta_Z1 = TofHit[0].beta;
- Beta_Z2 = TofHit[1].beta;
- Gamma1 = 1. / sqrt(1 - Beta_Z1 * Beta_Z1);
- Gamma2 = 1. / sqrt(1 - Beta_Z2 * Beta_Z2);
-
- AoQ1 = Brho1 / (3.10761 * Beta_Z1 * Gamma1);
- AoQ2 = Brho2 / (3.10761 * Beta_Z2 * Gamma2);
-
- A1 = AoQ1 * iZ1;
- A2 = AoQ2 * iZ2;
-
- // *** Filling the Fission Fragment Tree *** //
- SofFF->SetTOF(TofHit[0].tof);
- SofFF->SetTOF(TofHit[1].tof);
- SofFF->SetTofPosX(TofHit[0].x);
- SofFF->SetTofPosX(TofHit[1].x);
- SofFF->SetTofPosY(TofHit[0].y);
- SofFF->SetTofPosY(TofHit[1].y);
- SofFF->SetPlastic(TofHit[0].plastic);
- SofFF->SetPlastic(TofHit[1].plastic);
-
- SofFF->SetPosX1(TofHit[0].x2twim);
- SofFF->SetPosX1(TofHit[1].x2twim);
- SofFF->SetPosX2(TofHit[0].x2);
- SofFF->SetPosX2(TofHit[1].x2);
- SofFF->SetPosX3(TofHit[0].x3);
- SofFF->SetPosX3(TofHit[1].x3);
- SofFF->SetPosY3(TofHit[0].y3);
- SofFF->SetPosY3(TofHit[1].y3);
-
- SofFF->SetThetaIn(TofHit[0].theta_in/deg);
- SofFF->SetThetaIn(TofHit[1].theta_in/deg);
- SofFF->SetThetaOut(TofHit[0].theta_out/deg);
- SofFF->SetThetaOut(TofHit[1].theta_out/deg);
-
- SofFF->SetBeta(Beta_Z1);
- SofFF->SetBeta(Beta_Z2);
- SofFF->SetGamma(Gamma1);
- SofFF->SetGamma(Gamma2);
- SofFF->SetiZ(iZ1);
- SofFF->SetiZ(iZ2);
- SofFF->SetZ(Z1);
- SofFF->SetZ(Z2);
- SofFF->SetAoQ(AoQ1);
- SofFF->SetAoQ(AoQ2);
- SofFF->SetA(A1);
- SofFF->SetA(A2);
- SofFF->SetBrho(Brho1);
- SofFF->SetBrho(Brho2);
-
- SofFF->SetDT(TofHit[0].DT);
- SofFF->SetDT(TofHit[1].DT);
- SofFF->SetSection(TofHit[0].section);
- SofFF->SetSection(TofHit[1].section);
-
- SofFF->SetZsum(Zsum);
- SofFF->SetiZsum(iZsum);
-
- SofFF->SetFlightPath(TofHit[0].flight_path);
- SofFF->SetFlightPath(TofHit[1].flight_path);
+ SofFF->SetZsum(TofHit[0].Z+TofHit[1].Z);
+ SofFF->SetiZsum(TofHit[0].iZ+TofHit[1].iZ);
}
}
}
diff --git a/Projects/s455/RunToTreat.txt b/Projects/s455/RunToTreat.txt
index 96e1b4884bedc4fb8ffeaf46974e7f539bd1bf8e..aebd6fbe7df474485f907f82f41f1e3d87190294 100644
--- a/Projects/s455/RunToTreat.txt
+++ b/Projects/s455/RunToTreat.txt
@@ -6,6 +6,7 @@ RootFileName
%/media/sofia/s455/raw/run_raw_0368.root
%/media/sofia/s455/raw/run_raw_0369.root
%/media/sofia/s455/raw/run_raw_0419.root
+ %/media/sofia/s455/raw/run_raw_0421.root
%/media/sofia/s455/raw/run_raw_0422.root
%% 216Th