Truncate humongous eloss files, clean up Makevertex a little

- Changed Kinematics.h to use TLorentzVector - doesn't affect Ex calculations at all
	- 27Al, 17F analyses fleshed out to identify vertexz regions sensitive to a, d, p ejectiles
	- Better handling of dE using detector hyperboloids, dE4 introduced. Marginal improvement over dE3. This is a function inside PC_Stepladder_Correction.h which needs to be renamed.
	- Anode-wire #23 seems to have disproportionally large # of counts compared to the other wires.
	- Move towards wire-wise dE/E gates. Not yet fully implemented.
	- This push sees some Ex plots that appear to be from (a,p) on 27Al(a,p). 17F(a,p) is a lot less clear-cut.
	- Updated beam starting energies to improved guesses, for p+a, 17F and 27Al data.

Armory/ClassPW.h

@@ -30,7 +30,7 @@ struct PWHitInfo
 
 struct Coord
 {
-  float x, y, z;
+  double x, y, z;
   Coord() : x(0), y(0), z(0) {}
   Coord(const TVector3 &vec)
   {
@@ -260,7 +260,32 @@ inline void PW::PrintGeometry() {
 	}
 }
 
-inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair<TVector3, TVector3> awire, double sx3phi_radian) {
+inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair<TVector3, TVector3> awire, double phi) //Faster function
+{
+  const TVector3& a1 = awire.first;
+  const TVector3& a2 = awire.second;
+  const double s = TMath::Sin(phi), c = TMath::Cos(phi);
+  const double dx = a2.X() - a1.X(), dy = a2.Y() - a1.Y();
+  const double t = (a1.Y()*c - a1.X()*s) / (dx*s - dy*c); // denom != 0: wires always twist
+
+  auto nearerEndpoint = [&]() -> TVector3 {
+    auto dphi = [&](const TVector3& p) {
+      return TMath::Abs(TVector2::Phi_mpi_pi(phi - p.Phi()));
+    };
+    return dphi(a1) <= dphi(a2) ? a1 : a2;
+  };
+
+  if (t < 0.0 || t > 1.0)
+    return nearerEndpoint();
+
+  const TVector3 hit = a1 + t * (a2 - a1);
+  if (hit.X()*c + hit.Y()*s <= 0.0)  // wrong half-plane (anti-phi side)
+    return nearerEndpoint();
+
+  return hit;
+}
+
+/*inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair<TVector3, TVector3> awire, double sx3phi_radian) {
     // 1. Get wire geometry
     TVector3 a1 = awire.first;  // Top of the wire
     TVector3 a2 = awire.second; // Bottom of the wire
@@ -296,7 +321,7 @@ inline TVector3 PW::getClosestWirePosAtWirePhi(std::pair<TVector3, TVector3> awi
         }
     }
     return best_pcz_intersect;
-}
+}*/
 
 
 inline std::vector<std::vector<std::tuple<int,double,double>>>
@@ -365,6 +390,13 @@ PW::Make_Clusters(std::unordered_map<int,std::tuple<int,double,double>> wireEven
   } 	*/
 }
 
+
+/*
+ * inline TVector3 getRawCrossover(int anodeIndex, int cathodeIndex1, int cathodeIndex2, double cathodeE1, double cathodeE2) {
+ *		return (Crossover[anodeIndex][cathodeIndex1][0]*cathodeE1 + Crossover[anodeIndex][cathodeIndex2][0]*cathodeE2)/(cathodeE1+cathodeE2);
+ * }
+ * 
+ * */
 inline std::tuple<std::pair<TVector3, TVector3>, double, double, double>
 PW::GetPseudoWire(const std::vector<std::tuple<int,double,double>>& cluster, std::string type) {
 	std::pair<TVector3,TVector3> avgvec = std::pair(TVector3(0,0,0),TVector3(0,0,0));
@@ -417,6 +449,32 @@ PW::GetPseudoWire(const std::vector<std::tuple<int,double,double>>& cluster, std
 	return std::tuple(avgvec, sumEnergy, maxEnergy, tsMaxEnergy); 
 }
 
+#include <vector>
+#include <algorithm>
+#include <TGraph.h>
+
+double FastEval(TGraph* graph, double x) {
+    Double_t* x_arr = graph->GetX();
+    Double_t* y_arr = graph->GetY();
+    Int_t n = graph->GetN();
+
+    // Ensure array is sorted before using binary search
+    auto it = std::lower_bound(x_arr, x_arr + n, x);
+
+    if (it == x_arr) return y_arr[0];
+    if (it == x_arr + n) return y_arr[n - 1];
+
+    int idx = std::distance(x_arr, it);
+    double x1 = x_arr[idx - 1];
+    double x2 = x_arr[idx];
+    double y1 = y_arr[idx - 1];
+    double y2 = y_arr[idx];
+
+    // Linear interpolation
+    return y1 + (x - x1) * (y2 - y1) / (x2 - x1);
+}
+
+
 inline std::tuple<TVector3,double,double,double,double,double,double,double> PW::FindCrossoverProperties(const std::vector<std::tuple<int,double,double>>& a_cluster,
 																				 const std::vector<std::tuple<int,double,double>>& c_cluster) {
 	//std::pair<TVector3, TVector3> apwire = GetPseudoWire(a_cluster,"ANODE",anodeSumE);

Armory/HistPlotter.h

@@ -15,7 +15,6 @@
 #include <utility>
 #include <fstream>
 #include <sstream>
-#include <map>
 #include <unordered_map>
 #include <set>
 #include <TGraphErrors.h>
@@ -24,7 +23,7 @@ class HistPlotter {
 private:
     long long barrier_count, barrier_limit; //meant to keep track of how often to call FillN() on histograms
     enum {TFILE, TMEMFILE} filetype;
-    std::map<std::string,TObject*> oMap; //!< Maps std::string to all TH1, TH2 objects in the class
+    std::unordered_map<std::string,TObject*> oMap; //!< Maps std::string to all TH1, TH2 objects in the class
     std::unordered_map<std::string,TObject*> cutsMap; //!< Maps std::string to TCutG objects held by the class
     std::set<std::string> folderList; //!< List of all folder names used to nest objects
     std::unordered_map<TObject*,std::string> foldersForObjects; //!< Map that returns the folder corresponding to the object whose pointer is specified
@@ -32,8 +31,8 @@ private:
     TMemFile *omfile=nullptr; //!< TFile pointer for the output memfile
 
     //Caches to permit FillN() calls
-    std::map<std::string, std::vector<double>> onedimcache;
-    std::map<std::string, std::pair<std::vector<double>, std::vector<double>>> twodimcache;
+    std::unordered_map<std::string, std::vector<double>> onedimcache;
+    std::unordered_map<std::string, std::pair<std::vector<double>, std::vector<double>>> twodimcache;
     inline void FillN_All_Histograms();
 public:
     HistPlotter(std::string outfile, std::string type);

Armory/Kinematics.h

@@ -143,10 +143,12 @@ double Kinematics::getExc(double t3, double angle3)
     theta4 = (180./M_PI)*TMath::ASin((p3/P4)*TMath::Sin(angle3*M_PI/180.));
 
     //recalculate everything other than angle with lowered Kinetic energy
-    P4 = TMath::Sqrt(T4*T4 + 2*m4*T4);
+    P4 = TMath::Sqrt(T4*T4 + 2*(m4+Q-Q0)*T4);
     gamma4 = T4/m4+1.;
     beta4 = TMath::Sqrt(1. - 1./(gamma4*gamma4));
     theta4 = (180./M_PI)*TMath::ASin((p3/P4)*TMath::Sin(angle3*M_PI/180.));
+    
+    //if(m2==m3) std::cout << "Q0 elastic: " << Q0 << " " << Q << "\n" ;
 
     return Q0 - Q;//Q0  = Q + Exc
 

Armory/PC_StepLadder_Correction.h

@@ -1,14 +1,4 @@
 #include <TF1.h>
-/*double model(double* x, double* p) {
-    double result = x[0];
-    double factor = 29.0;
-    double slope = 0.7;
-    if(TMath::Abs(x[0]) < 16.2) result=x[0]*slope;
-    else if(TMath::Abs(x[0]) < 49.8 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor;
-    else if(TMath::Abs(x[0]) < 85.2 ) result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*2;
-    else result=x[0]*slope+TMath::Sign(1.0,x[0])*factor*3;
-    return result;
-}*/
 
 double model_invert(double* y, double* p) {
     double result = y[0];
@@ -37,9 +27,46 @@ double model_a1c1(double* x, double* p) {
 
 double model_invert_a1c1(double *y, double *q) {
     double result=y[0];
-    return result+40;
+    return result;
 }
 
+inline std::tuple<TVector3,TVector3,double> find_PC_PathLength(const TVector3& x1, const TVector3& x2) {
+    /*
+        Function that finds the path length between anode and cathode surfaces, both one-sheet hyperboloids of form (x*x+y*y)/(a*a) - (z*z)/(c*c) = 1
+        * path length found for a given particle moving along a certain direction from x1 to x2
+        * Typical arguments here will be x1=r_rhoMin, x2=qqqevent.pos or sx3event.pos
+    */
+    TVector3 dx = x2-x1;// direction vector
+    double t2 = 1.0; //The value of 't' at the destination point, by definition: t=(z(t)-z0)/dz 
+    auto onesheet_hyperboloid_intersect = [&](double a, double c) {
+		auto A = pow(dx.Perp(),2)/(a*a) - pow(dx.Z(),2)/(c*c);
+		auto B = 2*(dx.X()*x1.X()+dx.Y()*x1.Y())/(a*a) - 2*(dx.Z()*x1.Z())/(c*c);
+		auto C = pow(x1.Perp(),2)/(a*a) - pow(x1.Z(),2)/(c*c) - 1.0;
+		double disc = B*B - 4*A*C;
+		if(disc<0) 
+			return TVector3(0,0,54321);
+		else {
+			double tsol1 = (-B + TMath::Sqrt(disc))/(2*A);
+			double tsol2 = (-B - TMath::Sqrt(disc))/(2*A);
+			if(tsol1 >= 0  && tsol1 <= t2) 
+				return x1+tsol1*dx;
+			else if(tsol2>=0 && tsol2 <= t2)
+				return x1+tsol2*dx;
+			else
+				return TVector3(0,0,54321);
+		}
+	};
+
+	//TODO: Magic numbers here describing waist 'a', and flare 'c' will need better treatment. 
+	//Currently, these are derived by fitting the crossover points to R^2/a^2 - z^2/c^2 = 1 for anodes
+	// Cathode a, c values are found by scaling up the anode waist by 43/37, the ratio of the outermost radii
+	TVector3 anode_intersect = onesheet_hyperboloid_intersect(32.0429,301.895);
+	TVector3 cathode_intersect = onesheet_hyperboloid_intersect(37.239045,301.895);
+	if(anode_intersect.Z()!=54321 && cathode_intersect.Z()!=54321)
+		return std::tuple(cathode_intersect,anode_intersect,(cathode_intersect-anode_intersect).Mag()*0.1);
+	else
+		return std::tuple(TVector3(0,0,0), TVector3(0,0,0), 54321);
+}
 
 /*void testmodel() {
     TF1 eqline("x","x",-200,200);

results/compare.C

@@ -1,5 +1,5 @@
 {
-    TGraph kinematics("a(p,p)a_kinematics_7MeV_p.txt","%*lf %lf %*lf %lf");
+    TGraph kinematics("a(p,p)a_kinematics_6.88MeV_p.txt","%*lf %lf %*lf %lf");
     kinematics.Scale(4.0);
     TGraph kin2(kinematics.GetN(), kinematics.GetY(),kinematics.GetX());
     //TFile fin("4He(pp)_candidate_kinematic_curve.root");

results/compareE.C

@@ -1,5 +1,5 @@
 {
-    TGraph kinematics("a(p,p)a_kinematics_7MeV_p.txt","%*lf %*lf %lf %lf");
+    TGraph kinematics("a(p,p)a_kinematics_6.88MeV_p.txt","%*lf %*lf %lf %lf");
     kinematics.Scale(4.0);
     kinematics.SetLineColor(kBlue);
     kinematics.SetLineWidth(2.);

results/compareE_pcalpha.C

@@ -1,7 +1,7 @@
 {
     TCanvas c;
     //c.SetLogy(kTRUE);
-    TGraph kinematics("a(p,p)a_kinematics_7MeV_p.txt","%*lf %*lf %lf %lf");
+    TGraph kinematics("a(p,p)a_kinematics_6.88MeV_p.txt","%*lf %*lf %lf %lf");
     kinematics.Scale(4.0);
     kinematics.SetLineColor(kBlue);
     kinematics.SetLineWidth(2.);

results/haddnow

@@ -1,6 +1,6 @@
 rm out.root
-#hadd out.root ../results_run15.root ../results_run18.root ../results_run19.root ../results_run17.root ../results_run20.root ../results_run21.root ../results_run22.root
+hadd -k out.root ../results_run15.root ../results_run18.root ../results_run19.root ../results_run17.root ../results_run20.root ../results_run21.root ../results_run22.root
 #hadd out.root ../results_run015.root ../results_run018.root ../results_run019.root ../results_run017.root ../results_run020.root ../results_run021.root ../results_run022.root
-hadd out.root ../results_run018.root ../results_run019.root ../results_run017.root ../results_run020.root ../results_run021.root ../results_run022.root
+#hadd out.root ../results_run018.root ../results_run019.root ../results_run017.root ../results_run020.root ../results_run021.root ../results_run022.root
 
 #hadd out.root  ../results_run18.root ../results_run19.root ../results_run17.root ../results_run20.root ../results_run21.root ../results_run22.root

run_17F.sh

@@ -1,6 +1,10 @@
 rm results_run*.root
 export DATASET="17F"
+export flip180="0"
+export flipa=0
 export reactiondata=1
+export anode_offset=0
+export pressure_in_torr=350
 root -l -q -x -e ".L MakeVertex.C++g"
 rm 17F_output/*.root
 
@@ -19,15 +23,16 @@ function run_once() {
     file_exists=$(test -f /home/sud/Desktop/Software2/ANASEN_analysis/data/17F_fsu_files/Run_"$wrun"_mapped.root)
     if [[ $file_exists -ne 0 ]]; then return; fi
     root -q -l -b -x /home/sud/Desktop/Software2/ANASEN_analysis/data/17F_fsu_files/Run_"$wrun"_mapped.root -e 'gInterpreter->ProcessLine("#define ADD_NEW_BRANCHES 1");' -e $(printf 'tree->Process("MakeVertex.C+","analyzed_run%s.root")' "$wrun");
-    mv analyzed_run$wrun.root results_run$wrun.root;
+    mv analyzed_run$wrun.root 17F_output/results_run$wrun.root;
 }
 
 export -f run_once
 #run_once 351
-#parallel -j 6 --ctag run_once {1} ::: {350..400}
-parallel -j 6 --ctag run_once {1} ::: {350..360}
+time parallel -j 6 --ctag run_once {1} ::: {350..400}
+#parallel -j 6 --ctag run_once {1} ::: {380..400}
+#parallel -j 6 --ctag run_once {1} ::: {350..360}
 rm output_17F.root
-hadd -j 4 -k output_17F.root results_run3*.root
+hadd -j 4 -k output_17F.root 17F_output/results_run3*.root
 
 unset souce_vertex
 unset DATASET

run_27Al.sh

@@ -4,6 +4,7 @@ export flip180="0"
 export flipa=0
 export reactiondata=1
 export anode_offset=0
+export pressure_in_torr=250
 root -l -q -x -e ".L MakeVertex.C++"
 rm 27Al_output/*.root
 #declare -i run=28
@@ -29,14 +30,19 @@ function run_once() {
 }
 
 export -f run_once
-time parallel -j 2 --ctag run_once {1} ::: {50..59}
-time parallel -j 2 --ctag run_once {1} ::: 62 63 66 67 68 73 74
-#time parallel -j 2 --ctag run_once {1} ::: {78..89}
-
-
+function run_all() {
+    time parallel -j 6 --ctag run_once {1} ::: {50..59}
+    time parallel -j 6 --ctag run_once {1} ::: 62 63 66 67 73 74
+    time parallel -j 1 --ctag run_once {1} ::: 68
+    time parallel -j 6 --ctag run_once {1} ::: {78..89}
+}
+export -f run_all
+time run_all
+rm output.root
 hadd -k -j 4 output.root 27Al_output/results_run*.root
 mv output.root output_27Al.root
 rootbrowse output_27Al.root
 unset souce_vertex
 unset DATASET
 unset flip180
+unset pressure_in_torr

run_sx3.sh

@@ -1,5 +1,5 @@
 #Alpha runs at different spacer positions
-#rm results_run*.root
+rm results_run*.root
 export anode_offset=0
 export cathode_offset=0
 export DATASET="27Al"
@@ -14,6 +14,7 @@ root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_005_mapped.root -e 'tree-
 root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_006_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run06.root;
 root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_007_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run07.root;
 root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_008_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run08.root;
+exit
 fi
 
 #exit
@@ -34,12 +35,11 @@ unset timecut_low
 fi
 
 #protons+gas, 27Al
-if [[ 1 -eq 0 ]] ; then
+if [[ 1 -eq 1 ]] ; then
 #export anode_offset=0
-#export source_vertex=-200.0; #put the 'source' on the entrance window
 #export source_vertex=-57.28; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_015_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run15.root;
 export source_vertex=-135.68; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_017_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run17.root;
-exit
+#exit
 root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_018_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run18.root;
 root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_019_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run19.root;
 root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_020_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run20.root;
@@ -48,7 +48,7 @@ root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_022_mapped.root -e 'tree->Pr
 exit
 fi
 
-if [[ 1 -eq 1 ]]; then
+if [[ 1 -eq 0 ]]; then
 function run_once() {
     wrun=$(printf "%03d" $1)
     file_exists=$(test -f ../ANASEN_analysis/data/27Al_Data/Run_"$wrun"_mapped.root)
@@ -102,7 +102,7 @@ exit
 fi
 
 #17F alpha run with gas
-if [[ 1 -eq 0 ]]; then
+if [[ 1 -eq 1 ]]; then
 export source_vertex=53.44; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run18.root;
 export source_vertex=14.24; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_019_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run19.root;
 export source_vertex=-24.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run20.root;