From 5c0b3000787901be5043b136ae76a6e38adf3640 Mon Sep 17 00:00:00 2001
From: Charlie Paxman <cp00474@surrey.ac.uk>
Date: Tue, 8 Nov 2022 13:39:13 +0000
Subject: [PATCH] *e793s

---
 Projects/e793s/Analysis.cxx            |  17 ++-
 Projects/e793s/Analysis.h              |   2 +-
 Projects/e793s/macro/CS2_dt.h          |  33 +++---
 Projects/e793s/macro/DrawPlots.h       | 150 +++++++++++++++++++++++++
 Projects/e793s/macro/KnownPeakFitter.h |  15 ++-
 5 files changed, 194 insertions(+), 23 deletions(-)

diff --git a/Projects/e793s/Analysis.cxx b/Projects/e793s/Analysis.cxx
index c04cbfdd9..aa1eabc63 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 071c320aa..150d05057 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 a479ea275..87614b9b6 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 4456c7e58..44eac5740 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 42b6a3ccd..f4f66e89c 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);
+}
+
-- 
GitLab