From c7022178afc73df28e2c3c1ab41b5db24d59bbc3 Mon Sep 17 00:00:00 2001
From: Charlie Paxman <cp00474@surrey.ac.uk>
Date: Fri, 27 May 2022 14:53:45 +0100
Subject: [PATCH] * e793s - add eloss to diagnostic output

---
 .../e793s/macro/BeamSpot/MinimizeBeamSpot.cxx | 33 ++++++++----
 .../e793s/macro/BeamSpot/MinimizeBeamSpot.h   | 53 +++++++++++++++----
 2 files changed, 66 insertions(+), 20 deletions(-)

diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
index 8c601ccca..35521abf0 100755
--- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.cxx
@@ -22,6 +22,7 @@ double devE(const double * parameter) {
   h5 -> Reset();
   h7 -> Reset();
   hT -> Reset();
+  hEL -> Reset();
 
   //Now that initial range is wide, crop to single peak
   if(refE==0.143){
@@ -75,13 +76,14 @@ double devE(const double * parameter) {
       0.5 * parameter[3] * micrometer, //pass through half target
       ThetaTarget                      //angle leaving target
     );
-
+    
 //cout << Energy << endl;
     //Final value of Ex
     double Ex = reaction.ReconstructRelativistic(Energy, ThetaTarget);
 
     //Fill histograms with
 //if(ThetaTarget/deg<130.){ // TESTING
+//if(ThetaTarget/deg>130.){ // TESTING
     if(allButMG3){
       if(detnum[i]!=3){
         h -> Fill(Ex);
@@ -93,6 +95,7 @@ double devE(const double * parameter) {
 //}
     DetectorSwitch(detnum[i], Ex);
     hT -> Fill(ThetaTarget/deg,Ex);
+    hEL -> Fill(ThetaTarget/deg,Energy-energy[i]);//ELoss in target and Al
   }
 
   //Initilise, Draw & Fit histograms
@@ -103,11 +106,11 @@ double devE(const double * parameter) {
 
   /*** Minimize by one peak ***/
 /**/
-  double multiplier = 0.10; //0.08;
-  //double metric = abs(FitResultMatrix[mgSelect][0]-refE) + abs(multiplier*FitResultMatrix[mgSelect][2]);
+  double multiplier = 0.05; //0.08;
+  double metric = abs(FitResultMatrix[mgSelect][0]-refE) + abs(multiplier*FitResultMatrix[mgSelect][2]);
   //double metric = abs(FitResultMatrix[mgSelect][0]-0.143) + abs(multiplier*FitResultMatrix[mgSelect][2]) + abs(FitResultMatrix[mgSelect][5]-1.410) + abs(multiplier*FitResultMatrix[mgSelect][6]);
   //double metric = abs(FitResultMatrix[mgSelect][0]-0.143) + abs(multiplier*FitResultMatrix[mgSelect][2]) + abs(FitResultMatrix[mgSelect][5]-1.410) + abs(multiplier*FitResultMatrix[mgSelect][6])+ abs(FitResultMatrix[mgSelect][7]-1.980) + abs(multiplier*FitResultMatrix[mgSelect][8]) ;
-  double metric = abs(FitResultMatrix[mgSelect][0]-0.143) + abs(multiplier*FitResultMatrix[mgSelect][2]);
+  //double metric = abs(FitResultMatrix[mgSelect][7]-1.981) + abs(multiplier*FitResultMatrix[mgSelect][8]); ;
 
 /**/
 
@@ -151,9 +154,8 @@ double devE(const double * parameter) {
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 void MinimizeBeamSpot() {
 
-  filename = "XYZE_GateOn1838_11Apr22.txt"; refE = 1.981; // the energy of the selected states
-  
-  //filename = "XYZE_GateOn0143_11Apr22.txt"; refE = 0.143; // the energy of the selected states
+  //filename = "XYZE_GateOn1838_11Apr22.txt"; refE = 1.981; // the energy of the selected states
+  filename = "XYZE_GateOn0143_11Apr22.txt"; refE = 0.143; // the energy of the selected states
 
 
   //Read data in
@@ -214,7 +216,16 @@ void MinimizeBeamSpot() {
     //-5.088554, +1.476001, +0.342552, +2.586484 
     //-5.088554, +1.476001, +0.236014, +2.523850 
     //-3.989139, 0.99116, 0.115426, 2.586484 
-    -3.989139, +0.991160, +0.587660, +2.586484
+    //-3.989139, +0.991160, +0.587660, +2.586484
+
+    //-3.989139, +0.991160, +0.587660, +2.586484
+    //-3.9164, +0.0550, 0.0, 2.6
+    //-3.9164, +0.0550, 0.0, +2.838327
+    //-3.912110, +0.083580, +0.000000, +2.838327
+    //-3.912110, +0.083580, +0.000000, +2.6
+    //-3.9164, +0.0550, +0.0, 2.3
+    //-3.9164, +0.0550, +0.0, +2.591512
+    -5.334864, -0.698900, +0.0, +2.591512
 
   };
 
@@ -239,15 +250,15 @@ void MinimizeBeamSpot() {
   minim -> SetFunction(func);
 
   //Assign variable limits
-  //minim -> SetLimitedVariable(0, "X", parameter[0], 0.20, -8.0, -0.0);
+  //minim -> SetLimitedVariable(0, "X", parameter[0], 0.10, -8.0, -0.0);
   minim -> SetFixedVariable(0, "X", parameter[0]); 
-  //minim -> SetLimitedVariable(1, "Y", parameter[1], 0.20, -6.0, +6.0);
+  //minim -> SetLimitedVariable(1, "Y", parameter[1], 0.10, -6.0, +6.0);
   minim -> SetFixedVariable(1, "Y", parameter[1]);
   //minim -> SetLimitedVariable(2, "Z", parameter[2], 0.05, 0.90, +1.10);
   //minim -> SetLimitedVariable(2, "Z", parameter[2], 0.10, -2.00, +2.00);
   minim -> SetFixedVariable(2, "Z", parameter[2]);
   //minim -> SetLimitedVariable(3, "T", parameter[3], 0.05, +0.9, +1.5); // ELASTICS, 1.2(3)
-  //minim -> SetLimitedVariable(3, "T", parameter[3], 0.20, +2.05, +3.25); // ELASTICS, 2.65(59)
+  //minim -> SetLimitedVariable(3, "T", parameter[3], 0.05, +2.05, +3.25); // ELASTICS, 2.65(59)
   minim -> SetFixedVariable(3, "T", parameter[3]);
 
 
diff --git a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
index 4f772a902..4fef19679 100755
--- a/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
+++ b/Projects/e793s/macro/BeamSpot/MinimizeBeamSpot.h
@@ -32,7 +32,7 @@ int writeCount = 0;
 //Histograms
 string filename;
 double refE; // the energy of the selected states
-static auto h = new TH1D("h","All MG#'s", 400,-1.0,3.0);
+static auto h = new TH1D("h","All MG#'s", 200,-1.0,3.0);
 static auto h1 = new TH1D("h1","Individual MG#'s", 80,-1.0,3.0);
 static auto h2 = new TH1D("h2","h2", 80,-1.0,3.0);
 static auto h3 = new TH1D("h3","h3", 80,-1.0,3.0);
@@ -40,6 +40,10 @@ static auto h4 = new TH1D("h4","h4", 80,-1.0,3.0);
 static auto h5 = new TH1D("h5","h5", 80,-1.0,3.0);
 static auto h7 = new TH1D("h7","h7", 80,-1.0,3.0);
 static auto hT = new TH2F("hT","hT", 60,100.,160.,80,-1.0,3.0);
+static auto hEL = new TH2F("hEL","hEL", 60,100.,160.,2000,0.0,0.1);
+
+
+
 
 Double_t f_full(Double_t *x, Double_t *par) {
   float xx = x[0];
@@ -55,6 +59,13 @@ Double_t f_full(Double_t *x, Double_t *par) {
   return result;
 }
 
+Double_t f_one(Double_t *x, Double_t *par) {
+  float xx = x[0];
+  double result, norm;
+  result = (par[3]/(par[1]*sqrt(2*pi)))
+	      * exp(-0.5*pow((xx-par[2])/par[1],2));
+  return result;
+}
 
 //static auto hT = new TH2F("hT","hT", 20,100.,160.,20,-1.0,1.0);
 ////////////////////////////////////////////////////////////////////////////////
@@ -229,8 +240,27 @@ void DrawOneHistogram(TH1D* hist, int mg, int colour, int fill, double *FitResul
   } else {
     settings = "RBWQSN";
   }  
-
-
+/*
+  //TFitResultPtr fit;
+  //if(abs(refE-0.143)<0.05){
+    TF1 *full = new TF1("fitThreePeaks", f_full, -1.0, +2.6, (int) 1+(3*3));
+    for(int i=0; i<3; i++) {
+      full->SetParameter((i*3)+1,0.14);
+      full->SetParameter((i*3)+3,1e3);
+    }
+    full->SetParameter((0*3)+2,0.143);
+    full->SetParameter((1*3)+2,1.410);
+    full->SetParameter((2*3)+2,1.981);
+    //fit = hist->Fit(full, settings, "", -1.0, +2.6);
+    TFitResultPtr fit = hist->Fit(full, settings, "", -1.0, +2.6);
+  //} else {
+  //  TF1 *one = new TF1("fitOne", f_one, -1.0, +2.6, (int) 1+(3*3));
+  //  one->SetParameter(1,0.14);
+  //  one->SetParameter(3,1e3);
+  //  one->SetParameter(2,refE);
+  //  fit = hist->Fit(one, settings, "", -1.0, +2.6);
+  //}
+*/
   TF1 *full = new TF1("fitThreePeaks", f_full, -1.0, +2.6, (int) 1+(3*3));
   for(int i=0; i<3; i++) {
     full->SetParameter((i*3)+1,0.14);
@@ -241,10 +271,6 @@ void DrawOneHistogram(TH1D* hist, int mg, int colour, int fill, double *FitResul
   full->SetParameter((2*3)+2,1.981);
 
   TFitResultPtr fit = hist->Fit(full, settings, "", -1.0, +2.6);
-
-  
-
-  //TFitResultPtr fit = hist->Fit("gaus",settings); //N=stop drawing, Q=stop writing
   FillMatrix(FitResultMatrixMG,fit);
 } 
 ////////////////////////////////////////////////////////////////////////////////
@@ -254,7 +280,8 @@ void InitiliseCanvas(double FitResultMatrix[7][9]){
   TCanvas *canv = new TCanvas("canv","Ex Histograms",20,20,1600,800);
   gStyle->SetOptStat(0);
   //canv->Divide(2,1,0.005,0.005,0);
-  canv->Divide(3,1,0.005,0.005,0);
+  //canv->Divide(3,1,0.005,0.005,0);
+  canv->Divide(2,2,0.005,0.005,0);
   canv->cd(1)->SetLeftMargin(0.15);
   canv->cd(1)->SetBottomMargin(0.15);
   gPad->SetTickx();
@@ -263,11 +290,14 @@ void InitiliseCanvas(double FitResultMatrix[7][9]){
   canv->cd(2)->SetBottomMargin(0.15);
   gPad->SetTickx();
   gPad->SetTicky();
-
   canv->cd(3)->SetLeftMargin(0.15);
   canv->cd(3)->SetBottomMargin(0.15);
   gPad->SetTickx();
   gPad->SetTicky();
+  canv->cd(4)->SetLeftMargin(0.15);
+  canv->cd(4)->SetBottomMargin(0.15);
+  gPad->SetTickx();
+  gPad->SetTicky();
 
   //Histogram setup - Individual
   canv->cd(1);
@@ -329,6 +359,11 @@ void InitiliseCanvas(double FitResultMatrix[7][9]){
   l0143->SetLineColor(kRed);
   l0143->Draw("same");
 
+  canv->cd(4);
+  hEL->GetXaxis()->SetTitle("Theta (degrees)");
+  hEL->GetYaxis()->SetTitle("Energy loss in Target + Al [MeV]");
+  hEL->Draw();
+
   //Refresh
   gPad->Update();
 
-- 
GitLab