angular straggling

Armory/anasenMS.cpp

@@ -51,6 +51,13 @@ TGraph* LoadSigmaXVsEnergy(const std::string& filename) {
   return g;
 }
 
+TGraph* LoadSigmaAVsEnergy(const std::string& filename) { //load column 2 vs column 5 (Sigma_a_cm) from an E_vs_a_*.dat table
+  TGraph* g = new TGraph(filename.c_str(), "%*lg %lg %*lg %*lg %lg");
+  g->Sort(); // TGraph::Eval requires ascending x (Energy_MeV)
+  return g;
+}
+
+
 bool IsDeadAnode(int id){
     static std::set<int> dead = {}; // add dead anode IDs here, 0-23
     return dead.count(id);
@@ -198,6 +205,7 @@ int main(int argc, char **argv){
   
   TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = path length (cm), y = beam energy (MeV)
   TGraph* sigmaXBeam = LoadSigmaXVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_x (cm)
+  TGraph* sigmaABeam = LoadSigmaAVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_a (cm)
   // Build a temporary inverse (energy -> path) to locate the path at beamE.
   TGraph* elossBeamInverseRaw = new TGraph(elossBeam->GetN());
   for( int p = 0; p < elossBeam->GetN(); p++ ){
@@ -422,9 +430,10 @@ int main(int argc, char **argv){
   tree1->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
   tree1->Branch("rePhi1",     &rePhi1, "reconstucted_phi1/D");
 
-  double hitTheta, originalEnergy;
+  double hitTheta, originalEnergy, sigma_a;
   tree1->Branch("hitTheta", &hitTheta, "hitTheta/D");
   tree1->Branch("originalEnergy", &originalEnergy, "originalEnergy/D");
+  tree1->Branch("sigma_a", &sigma_a, "sigma_a/D");
 
   // reconstructed vertex Z from PW fit
   double z0;
@@ -471,6 +480,8 @@ int main(int argc, char **argv){
     detX = TMath::QuietNaN();
     detY = TMath::QuietNaN();
     detZ = TMath::QuietNaN();
+    hitTheta = TMath::QuietNaN();
+    originalEnergy = TMath::QuietNaN();
     sx3ID = -1;
     sx3Up = -1;
     sx3Dn = -1;
@@ -496,13 +507,11 @@ int main(int argc, char **argv){
     
     // transverse sampling range from the beam's distance straggle at this energy
     const double sigmaX_mm = std::max(0.0, sigmaXBeam->Eval(beamEnergy)) * 10.0; // cm -> mm
-    //vertexX = 2.0 * sigmaX_mm * gRandom->Rndm() - sigmaX_mm;
-    //vertexY = 2.0 * sigmaX_mm * gRandom->Rndm() - sigmaX_mm;
     double vertexRangeX = std::sqrt((vertexXRange[1] * vertexXRange[1]) + (sigmaX_mm * sigmaX_mm));
     double vertexRangeY = std::sqrt((vertexYRange[1] * vertexYRange[1]) + (sigmaX_mm * sigmaX_mm));
     //std::cout << "vertexRangeX: " << vertexRangeX << ", vertexRangeY: " << vertexRangeY << std::endl;
-    vertexX = gRandom->Gaus(0, vertexRangeX); // mean and standard deviation
-    vertexY = gRandom->Gaus(0, vertexRangeY); // mean and standard deviation
+    vertexX = gRandom->Gaus(0, vertexRangeX) + 5; // mean and standard deviation
+    vertexY = gRandom->Gaus(0, vertexRangeY) - 5; // mean and standard deviation
 
     //vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() +  vertexZRange[0];
 
@@ -518,8 +527,8 @@ int main(int argc, char **argv){
     //double beamEnergyLoss = elossBeam->Eval(0.0) - beamEnergy;
     //KEA = beamEnergy / beamA;
     //KEA = gRandom->Uniform(0, beamE);
-
-    transfer.SetIncidentEnergyAngle(KEA, 0, 0);
+    sigma_a = std::max(0.0, sigmaABeam->Eval(beamEnergy));
+    transfer.SetIncidentEnergyAngle(KEA, sigma_a, 0); //arguments are (kinetic energy, polar angle, azimuthal angle) of the incident particle in the lab frame
     transfer.CalReactionConstant();
 
     // isotropic CM direction
@@ -582,7 +591,7 @@ int main(int argc, char **argv){
       sx3Up    = sx3->GetChUp();
       sx3Dn    = sx3->GetChDn();
       sx3Bk    = sx3->GetChBk();
-      if(IsDeadSX3ChannelCombo(sx3ID, sx3Up, sx3Dn, sx3Bk)) sx3Up = -1, sx3Dn = -1, sx3Bk = -1; // mark as no hit if any SX3 channel is dead
+      if(IsDeadSX3ChannelCombo(sx3ID, sx3Up, sx3Dn, sx3Bk)) sx3Up = -1, sx3Dn = -1, sx3Bk = -1, sx3ID = -1; // mark as no hit if any SX3 channel is dead
       sx3ZFrac = sx3->GetZFrac();
       
       // apply intrinsic detector resolution to true SX3 hit position

Armory/testMacro.C

@@ -24,13 +24,13 @@
 		100, -1, 0
 	);
 	tree1->Draw("Ex:vZ>>hExVsVZ", "sx3ID >= 0", "box same");
-
+/*
 	new TCanvas();
 	tree1->Draw("originalEnergy:Tb", "sx3ID >=0", "");
 	TF1 *f_diagonal = new TF1("f_diagonal", "x", 0, 200); // assuming the range for Tb is 0 to 200
 	f_diagonal->SetLineColor(kRed);
 	f_diagonal->SetLineStyle(2);
-	f_diagonal->Draw("same");
+	f_diagonal->Draw("same");*/
 
 
 }

ELoss/EnergyLoss.py

@@ -171,18 +171,21 @@ def make_E_vs_x(
     e_u_init = emax_mev / mass_u
     sigma_E = np.zeros_like(E)
     sigma_x = np.zeros_like(E)
+    sigma_a = np.zeros_like(E)
     for i, depth_cm in enumerate(x):
         gas.density(rho_g_cm3).thickness(depth_cm * rho_g_cm3)
         projectile.T(e_u_init)
         result = catima.calculate(projectile, gas)
         sigma_E[i] = result.sigma_E
         sigma_x[i] = result.sigma_x
+        sigma_a[i] = result.sigma_a
 
     df = pd.DataFrame({
         "Distance_cm": x,
         "Energy_MeV": E,
         "Sigma_E_MeV": sigma_E,
-        "Sigma_x_cm": sigma_x
+        "Sigma_x_cm": sigma_x,
+        "Sigma_a_cm": sigma_a
     })
 
     outfile = get_loss_table_path(medium, label)