From d18c626107eebde450a0cbb76dbc9fc4de14da08 Mon Sep 17 00:00:00 2001
From: Charlie Paxman <cp00474@surrey.ac.uk>
Date: Fri, 28 May 2021 18:22:41 +0100
Subject: [PATCH] * Correcting angle for target energy loss Was previously set
 to assume that the particle leaves the target at 0degrees, therefore
 underestimated the energy loss in the target Now uses proper angle
 calculation

---
 Projects/e793s/Analysis.cxx              | 25 +++++++--------
 Projects/e793s/macro/BeamSpot/BeamSpot.C | 41 ++++++++++++++++--------
 2 files changed, 38 insertions(+), 28 deletions(-)

diff --git a/Projects/e793s/Analysis.cxx b/Projects/e793s/Analysis.cxx
index 630354c87..60fa8381b 100755
--- a/Projects/e793s/Analysis.cxx
+++ b/Projects/e793s/Analysis.cxx
@@ -174,9 +174,6 @@ void Analysis::TreatEvent(){
     ThetaM2Surface = HitDirection.Angle(- M2->GetTelescopeNormal(countMust2) );
     ThetaNormalTarget = HitDirection.Angle( TVector3(XBeam,YBeam,1) ) ;
 
-   // cout<<"Must2 Znormal:"<<M2 -> GetTelescopeNormal(countMust2).Z()<<endl;
-  //  cout<<"Must2 telescope:"<<M2->TelescopeNumber[countMust2]<<endl;
-
     /************************************************/
     // Part 2 : Impact Energy
     Energy = 0;
@@ -195,9 +192,7 @@ void Analysis::TreatEvent(){
     else
       Energy = Si_E_M2;
 
-
-    RawEnergy.push_back(Energy); //CPx ADDITION
-
+    RawEnergy.push_back(Energy);
 
     // Evaluate energy using the thickness
     elab_tmp = LightAl.EvaluateInitialEnergy(Energy, 0.4*micrometer, ThetaM2Surface);
@@ -243,30 +238,32 @@ void Analysis::TreatEvent(){
     Y.push_back( MG -> GetPositionOfInteraction(countMugast).Y());
     Z.push_back( MG -> GetPositionOfInteraction(countMugast).Z());
 
-    //ThetaMGSurface = HitDirection.Angle( TVector3(0,0,1) ) ;
     ThetaMGSurface = HitDirection.Angle( MG -> GetTelescopeNormal(countMugast) );
     ThetaNormalTarget = HitDirection.Angle( TVector3(XBeam,YBeam,-1) ) ;
 
     // Part 2 : Impact Energy
     Energy = elab_tmp = 0;
     Energy = MG->GetEnergyDeposit(countMugast);
-
     RawEnergy.push_back(Energy);
 
-    elab_tmp = LightAl.EvaluateInitialEnergy(Energy, 0.4*micrometer, ThetaMGSurface);
-   
-    elab_tmp = LightTarget.EvaluateInitialEnergy(elab_tmp, 0.5*TargetThickness, 0.);
-    //elab_tmp   = LightTarget.EvaluateInitialEnergy( elab_tmp ,TargetThickness*0.5, ThetaNormalTarget);
+    elab_tmp = LightAl.EvaluateInitialEnergy(
+		    Energy,              //particle energy after Al
+		    0.4*micrometer,      //thickness of Al
+		    ThetaMGSurface);     //angle of impingement
+    elab_tmp = LightTarget.EvaluateInitialEnergy(
+		    elab_tmp,            //particle energy after leaving target
+		    TargetThickness*0.5, //distance passed through target
+		    ThetaNormalTarget);  //angle of exit from target
     ELab.push_back(elab_tmp);
 
     // Part 3 : Excitation Energy Calculation
     Ex.push_back(reaction.ReconstructRelativistic(elab_tmp,thetalab_tmp));
     Ecm.push_back(elab_tmp*(AHeavy+ALight)/(4*AHeavy*cos(thetalab_tmp)*cos(thetalab_tmp)));
+    
     // Part 4 : Theta CM Calculation
-    ThetaCM.push_back(reaction.EnergyLabToThetaCM(elab_tmp, thetalab_tmp)/deg);
     ThetaLab.push_back(thetalab_tmp/deg);
-     
     PhiLab.push_back(philab_tmp/deg);
+    ThetaCM.push_back(reaction.EnergyLabToThetaCM(elab_tmp, thetalab_tmp)/deg);
 
     if(sizeMG==1){
       MG_T = MG->DSSD_T[0];
diff --git a/Projects/e793s/macro/BeamSpot/BeamSpot.C b/Projects/e793s/macro/BeamSpot/BeamSpot.C
index c60835c10..7d9eaa2cc 100755
--- a/Projects/e793s/macro/BeamSpot/BeamSpot.C
+++ b/Projects/e793s/macro/BeamSpot/BeamSpot.C
@@ -29,21 +29,23 @@ void BeamSpot(){
   double ELab = 0.0, Ex = 0.0;
   vector <double> Xd, Yd, Zd;      //Vector of particle direction. Calculated as Xp-Xb, Yp-Yb...
   ifstream MugastDataFile;
+  double ThetaNormalTarget;
+
 
   gErrorIgnoreLevel = kWarning; // Suppress ".pdf created" lines
 
   /*** ITERATIVE GRID CONTROLS ***/
   /***** pos varied as offset ****/
-  /**/ double xmin = +0.000;   /**/
-  /**/ double xmax = +0.100;   /**/
-  /**/ unsigned int xdiv = 10; /**/
+  /**/ double xmin = +0.020;   /**/
+  /**/ double xmax = +0.080;   /**/
+  /**/ unsigned int xdiv = 12; /**/
   /**/                         /**/
-  /**/ double ymin = -0.100;   /**/
-  /**/ double ymax = +0.100;   /**/
+  /**/ double ymin = -0.050;   /**/
+  /**/ double ymax = +0.050;   /**/
   /**/ unsigned int ydiv = 10; /**/
   /**/                         /**/
-  /**/ double zmin = -0.100;   /**/
-  /**/ double zmax = +0.100;   /**/
+  /**/ double zmin = -0.050;   /**/
+  /**/ double zmax = +0.050;   /**/
   /**/ unsigned int zdiv =  2; /**/
   /**/                         /**/
   /***** thick varied as %ge *****/
@@ -58,9 +60,9 @@ void BeamSpot(){
   /*******************************/
 
   // File name controls
-  const char* XYZE_file = "XYZE_gammaGated_Run63.txt";
-  const char* outputMetric = "output_Run63_metrics.txt";
-  const char* outputHisto = "output_Run63_histograms.root";
+  const char* XYZE_file = "XYZE_gammaGated_Full_TestThetaNormalTarget.txt";
+  const char* outputMetric = "output_Run63_metrics_ThetaNormal.txt";
+  const char* outputHisto = "output_Run63_histograms_ThetaNormal.root";
 
   // Calculate size of iteratve steps
   double xstp = (xmax-xmin)/ ((double) xdiv);
@@ -153,7 +155,8 @@ void BeamSpot(){
   	 	                  Yp[i] - beamSpot.Y(),   //Yd
 			          Zp[i] - beamSpot.Z() }; //Zd
 	  
-	    tempTheta = ELab = Ex = 0.0;
+            ThetaNormalTarget = tempTheta = ELab = Ex = 0.0;
+
 
 	    switch(DetNum[i]){
 	      case 1:
@@ -179,13 +182,23 @@ void BeamSpot(){
 	        return; // Exit code
 	    }
 
+            // Change beam spot vector to inverse beam direction vector
+            beamSpot.SetZ(-1.0);
+            ThetaNormalTarget = particleDir.Angle(beamSpot);
+
 	    //micrometer defined in NPSystemOfUnits.h
-	    ELab = LightAl.EvaluateInitialEnergy(Ep[i], 0.4*micrometer, tempTheta); 
-	    ELab = LightTarget.EvaluateInitialEnergy(ELab, 0.5*TargetThickness, 0.);
+	    ELab = LightAl.EvaluateInitialEnergy(
+			    Ep[i],               //energy after Al 
+			    0.4*micrometer,      //thickness of Al
+			    tempTheta);          //angle of impingement
+	    ELab = LightTarget.EvaluateInitialEnergy(
+			    ELab,                //energy after leaving target 
+			    0.5*TargetThickness, //pass through half target
+			    ThetaNormalTarget);  //angle leaving target
 
             // Change beam spot vector to beam direction vector
 	    beamSpot.SetZ(1.0);
-            Ex = reaction.ReconstructRelativistic( ELab, particleDir.Angle(beamSpot) );
+            Ex = reaction.ReconstructRelativistic(ELab, particleDir.Angle(beamSpot));
 
 	    // Fill Ex histograms
 	    tempHist->Fill(Ex);
-- 
GitLab