diff --git a/Armory/ANASEN_model.C b/Armory/ANASEN_model.C
index 72bc262..8f1e952 100644
--- a/Armory/ANASEN_model.C
+++ b/Armory/ANASEN_model.C
@@ -119,3 +119,5 @@ void ANASEN_model(int anodeID1 = -1, int anodeID2 = -1, int cathodeID1 = -1, int
   geom->SetVisLevel(4);
   worldBox->Draw("ogle");
 }
+
+
diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h
index 9cce5ae..0fd64bc 100755
--- a/Armory/ClassPW.h
+++ b/Armory/ClassPW.h
@@ -91,7 +91,8 @@ public:
   void ConstructGeo();
   void FindWireID(TVector3 pos, TVector3 direction, bool verbose = false);
   void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
-  void CalTrack2(TVector3 sx3Pos, TVector3 anodeInt, bool verbose = false);
+  //void CalTrack2(TVector3 sx3Pos, TVector3 anodeInt, bool verbose = false);
+  void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose);
 
   void Print()
   {
@@ -447,21 +448,41 @@ inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbo
 }
 
 
-inline void PW::CalTrack2(TVector3 siPos, TVector3 anodeInt, bool verbose)
+inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose)
 {
 
-  double mx, my;
-  double z;
-  mx = siPos.X() / (siPos.X() - anodeInt.X());
-  my = siPos.Y() / (siPos.Y() - anodeInt.Y());
-  z=siPos.Z() + mx * (anodeInt.Z() - siPos.Z());
-  // if (mx == my)
-  {
-    trackVec=TVector3(0,0,z);
-  }
+  trackPos = sx3Pos;
+
+  double p1 = TMath::Abs(hitInfo.nearestDist.first + gRandom->Gaus(0, sigmaA));
+  double p2 = TMath::Abs(hitInfo.nextNearestDist.first + gRandom->Gaus(0, sigmaA));
+  double fracA = p1 / (p1 + p2);
+  short anodeID1 = hitInfo.nearestWire.first;
+  short anodeID2 = hitInfo.nextNearestWire.first;
+  TVector3 shiftA1 = (An[anodeID2].first - An[anodeID1].first) * fracA;
+  TVector3 shiftA2 = (An[anodeID2].second - An[anodeID1].second) * fracA;
+
+  double q1 = TMath::Abs(hitInfo.nearestDist.second + gRandom->Gaus(0, sigmaC));
+  double q2 = TMath::Abs(hitInfo.nextNearestDist.second + gRandom->Gaus(0, sigmaC));
+  double fracC = q1 / (q1 + q2);
+  short cathodeID1 = hitInfo.nearestWire.second;
+  short cathodeID2 = hitInfo.nextNearestWire.second;
+  TVector3 shiftC1 = (Ca[cathodeID2].first - Ca[cathodeID1].first) * fracC;
+  TVector3 shiftC2 = (Ca[cathodeID2].second - Ca[cathodeID1].second) * fracC;
+
+  TVector3 a1 = An[anodeID1].first + shiftA1;
+  TVector3 a2 = An[anodeID1].second + shiftA2;
+
+  TVector3 c1 = Ca[cathodeID1].first + shiftC1;
+  TVector3 c2 = Ca[cathodeID1].second + shiftC2;
+
+  TVector3 n1 = (a1 - a2).Cross((sx3Pos - a2)).Unit();
+  TVector3 n2 = (c1 - c2).Cross((sx3Pos - c2)).Unit();
+
+  // if the handiness of anode and cathode revered, it should be n2 cross n1
+  trackVec = (n2.Cross(n1)).Unit();
 
   if (verbose)
-    printf("X slope = %f and Y slope = %f \n", mx, my);
+    printf("Theta, Phi = %f, %f \n", trackVec.Theta() * TMath::RadToDeg(), trackVec.Phi() * TMath::RadToDeg());
 }
 
 /*inline TVector3 PW::CalTrack3(TVector3 siPos, TVector3 anodeInt, bool verbose)
diff --git a/Armory/Makefile b/Armory/Makefile
index 2ff7b6f..6693ef3 100644
--- a/Armory/Makefile
+++ b/Armory/Makefile
@@ -23,10 +23,11 @@ Mapper	:	Mapper.cpp  ../mapping.h ClassDet.h
 		@echo "--------- making Mapper"
 		$(CC) $(COPTS) -o Mapper Mapper.cpp $(ROOTLIBS)
 
-# AnasenMS : constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h anasenMS.cpp 
-# 		@echo "--------- making ANASEN Monte Carlo"
-# 		$(CC) $(COPTS) -o AnasenMS anasenMS.cpp $(ROOTLIBS)
+AnasenMS : constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h anasenMS.cpp 
+		@echo "--------- making ANASEN Monte Carlo"
+		$(CC) $(COPTS) -o AnasenMS anasenMS.cpp $(ROOTLIBS)
 
 EventBuilder : EventBuilder.cpp ClassData.h fsuReader.h Hit.h
 		@echo "--------- making EventBuilder"
 		$(CC) $(COPTS) -o EventBuilder EventBuilder.cpp $(ROOTLIBS) 
+
diff --git a/Armory/anasenMS_root b/Armory/anasenMS_root
new file mode 100755
index 0000000..cdb335e
Binary files /dev/null and b/Armory/anasenMS_root differ
diff --git a/Armory/anasenMS_root.cpp b/Armory/anasenMS_root.cpp
new file mode 100644
index 0000000..5a8ad39
--- /dev/null
+++ b/Armory/anasenMS_root.cpp
@@ -0,0 +1,309 @@
+// add includes at top
+#include <TApplication.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TEveManager.h>
+#include <TEvePointSet.h>
+#include <TEveGeoNode.h>
+#include <TGeoManager.h>
+#include <vector>
+#include <mutex>
+
+#include "TRandom.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TStyle.h"
+#include "TCanvas.h"
+#include "TBenchmark.h"
+
+#include "ClassTransfer.h"
+#include "ClassAnasen.h"
+
+// expose to ROOT
+int Run(int nEvents=1000, const char* outFile=nullptr){
+  // Ensure TEve exists (create after geometry has been built if possible)
+  if(!gEve) TEveManager::Create();
+
+  // if a geometry has already been loaded by ANASEN_model.C, make sure it
+  // shows up in the TEve scene.  TEveManager::Create() normally pulls in
+  // gGeoManager, but we do it explicitly to be safe.  We must wrap the
+  // top node/volume in a TEveGeoTopNode (not pass a raw TGeoVolume).
+  if(gGeoManager){
+    // create a TEve wrapper around the top node
+    TEveGeoTopNode *top = new TEveGeoTopNode(gGeoManager, gGeoManager->GetTopNode());
+    gEve->AddElement(top);
+  }
+
+  // Reaction
+  TransferReaction transfer;
+
+  transfer.SetA(24,12, 0);
+  transfer.SetIncidentEnergyAngle(10, 0, 0);
+  transfer.Seta( 4, 2);
+  transfer.Setb( 1, 1);
+
+  //TODO add alpha source
+
+  std::vector<float> ExAList = {0};
+  std::vector<float> ExList = {0, 1, 2};
+
+  double vertexXRange[2] = {   -5,    5}; // mm
+  double vertexYRange[2] = {   -5,    5}; 
+  double vertexZRange[2] = { -100,  100}; 
+
+  double sigmaSX3_W = -1; // mm, < 0 use mid-point
+  double sigmaSX3_L = 3; // mm, < 0 use mid-point
+  double sigmaPW_A  = 0; // from 0 to 1.
+  double sigmaPW_C  = 0; // from 0 to 1.
+
+  //###################################################
+
+  printf("------------ Vertex :\n");
+  printf("X : %7.2f - %7.2f mm\n", vertexXRange[0], vertexXRange[1]);
+  printf("Y : %7.2f - %7.2f mm\n", vertexYRange[0], vertexYRange[1]);
+  printf("Z : %7.2f - %7.2f mm\n", vertexZRange[0], vertexZRange[1]);
+  printf("------------ Uncertainty :\n");
+  printf(" SX3 horizontal : %.1f\n", sigmaSX3_W);
+  printf(" SX3   vertical : %.1f\n", sigmaSX3_L);
+  printf("          Anode : %.1f mm\n", sigmaPW_A);
+  printf("        Cathode : %.1f mm\n", sigmaPW_C);
+  printf("        num_eve : %d \n",nEvents);
+  transfer.CalReactionConstant();
+
+  int nExA = ExAList.size();
+  int nEx  = ExList.size();
+
+  ANASEN * anasen = new ANASEN();
+  SX3 * sx3 = anasen->GetSX3();    
+  PW * pw = anasen->GetPW(); 
+
+  TString saveFileName = "SimAnasen1.root";
+  printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
+  TFile * saveFile = new TFile(saveFileName, "recreate");
+  TTree * tree = new TTree("tree", "tree");
+
+  double KEA;
+  tree->Branch("beamKEA",     &KEA, "beamKEA/D");
+
+  double thetaCM, phiCM;
+  tree->Branch("thetaCM", &thetaCM, "thetaCM/D");
+  tree->Branch("phiCM",     &phiCM, "phiCM/D");
+
+  double thetab, phib, Tb;
+  double thetaB, phiB, TB;
+  tree->Branch("thetab", &thetab, "thetab/D");
+  tree->Branch("phib",     &phib, "phib/D");
+  tree->Branch("Tb",         &Tb, "Tb/D");
+  tree->Branch("thetaB", &thetaB, "thetaB/D");
+  tree->Branch("phiB",     &phiB, "phiB/D");
+  tree->Branch("TB",         &TB, "TB/D");
+
+  int ExAID;
+  double ExA;
+  tree->Branch("ExAID", &ExAID, "ExAID/I");
+  tree->Branch("ExA",     &ExA, "ExA/D");
+
+  int ExID;
+  double Ex;
+  tree->Branch("ExID", &ExID, "ExID/I");
+  tree->Branch("Ex",     &Ex, "Ex/D");
+
+  double vertexX, vertexY, vertexZ;
+  tree->Branch("vX",     &vertexX, "VertexX/D");
+  tree->Branch("vY",     &vertexY, "VertexY/D");
+  tree->Branch("vZ",     &vertexZ, "VertexZ/D");
+
+  double sx3X, sx3Y, sx3Z;
+  tree->Branch("sx3X",     &sx3X, "sx3X/D");
+  tree->Branch("sx3Y",     &sx3Y, "sx3Y/D");
+  tree->Branch("sx3Z",     &sx3Z, "sx3Z/D");
+
+  int anodeID[2], cathodeID[2];
+  tree->Branch("aID", anodeID, "anodeID/I");
+  tree->Branch("cID", cathodeID, "cathodeID/I");
+
+  double anodeDist[2], cathodeDist[2];
+  tree->Branch("aDist",   anodeDist, "anodeDist/D");
+  tree->Branch("cDist",  cathodeDist, "cathodeDist/D");
+
+  int sx3ID, sx3Up, sx3Dn, sx3Bk;
+  double sx3ZFrac;
+  tree->Branch("sx3ID",   &sx3ID, "sx3ID/I");
+  tree->Branch("sx3Up",   &sx3Up, "sx3Up/I");
+  tree->Branch("sx3Dn",   &sx3Dn, "sx3Dn/I");
+  tree->Branch("sx3Bk",   &sx3Bk, "sx3Bk/I");
+  tree->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
+
+  double reTheta, rePhi;
+  tree->Branch("reTheta", &reTheta, "reconstucted_theta/D");
+  tree->Branch("rePhi",     &rePhi, "reconstucted_phi/D");
+
+  double reTheta1, rePhi1;
+  tree->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
+  tree->Branch("rePhi1",     &rePhi1, "reconstucted_phi1/D");
+
+  double z0;
+  tree->Branch("z0",     &z0, "reconstucted_Z/D");
+
+
+  //========timer
+  TBenchmark clock;
+  bool shown ;   
+  clock.Reset();
+  clock.Start("timer");
+  shown = false;
+
+  // Create a point set to show hits
+  TEvePointSet *pts = new TEvePointSet("hits");
+  pts->SetMarkerStyle(20);
+  pts->SetMarkerColor(kRed);
+  gEve->AddElement(pts);
+
+  // Optionally open output file/tree
+  TFile *fout = nullptr;
+  TTree *tout = nullptr;
+  std::vector<double> vx, vy, vz;
+  if(outFile){
+    fout = TFile::Open(outFile,"RECREATE");
+    tout = new TTree("evt","events");
+    tout->Branch("x",&vx);
+    tout->Branch("y",&vy);
+    tout->Branch("z",&vz);
+  }
+
+  // Simulation loop (replace with your sim code that fills vx,vy,vz per event)
+
+  for( int i = 0; i < nEvents ; i++){
+
+    ExAID = gRandom->Integer(nExA);
+    ExA = ExAList[ExAID];
+    transfer.SetExA(ExA);
+
+    ExID = gRandom->Integer(nEx);
+    Ex = ExList[ExID];
+    transfer.SetExB(Ex);
+
+    transfer.CalReactionConstant();
+
+    thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ; 
+    phiCM   = (gRandom->Rndm() - 0.5) * TMath::TwoPi();
+
+    //==== Calculate reaction
+    TLorentzVector * output = transfer.Event(thetaCM, phiCM);
+    TLorentzVector Pb = output[2];
+    TLorentzVector PB = output[3];
+
+    thetab = Pb.Theta() * TMath::RadToDeg();
+    thetaB = PB.Theta() * TMath::RadToDeg();
+
+    Tb = Pb.E() - Pb.M();
+    TB = PB.E() - PB.M();
+
+    phib = Pb.Phi() * TMath::RadToDeg();
+    phiB = PB.Phi() * TMath::RadToDeg();
+
+    vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() +  vertexXRange[0];
+    vertexY = (vertexYRange[1]- vertexYRange[0])*gRandom->Rndm() +  vertexYRange[0];
+    vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() +  vertexZRange[0];
+
+    TVector3 vertex(vertexX, vertexY, vertexZ);
+
+    TVector3 dir(1, 0, 0);
+    dir.SetTheta(thetab * TMath::DegToRad());
+    dir.SetPhi(phib * TMath::DegToRad());
+
+   
+    pw->FindWireID(vertex, dir, false);
+    sx3->FindSX3Pos(vertex, dir, false);   
+
+    PWHitInfo hitInfo = pw->GetHitInfo();
+
+    anodeID[0]   = hitInfo.nearestWire.first;
+    cathodeID[0] = hitInfo.nearestWire.second;
+    anodeID[1]   = hitInfo.nextNearestWire.first;
+    cathodeID[1] = hitInfo.nextNearestWire.second;
+
+    anodeDist[0]   = hitInfo.nearestDist.first;
+    cathodeDist[0] = hitInfo.nearestDist.second;
+    anodeDist[1]   = hitInfo.nextNearestDist.first;
+    cathodeDist[1] = hitInfo.nextNearestDist.second;
+
+    sx3ID = sx3->GetID();
+    if( sx3ID >= 0 ){
+      sx3Up    = sx3->GetChUp();
+      sx3Dn    = sx3->GetChDn();
+      sx3Bk    = sx3->GetChBk();
+      sx3ZFrac = sx3->GetZFrac();
+
+      //Introduce uncertaity
+      // TVector3 hitPos = sx3->GetHitPos();
+      TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
+
+      sx3X = hitPos.X();
+      sx3Y = hitPos.Y();
+      sx3Z = hitPos.Z();
+      
+      pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
+      reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi = pw->GetTrackPhi() * TMath::RadToDeg();
+
+      pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false);
+      reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg();
+
+      z0 = pw->GetZ0();
+
+    }else{
+      sx3Up = -1;
+      sx3Dn = -1;
+      sx3Bk = -1;
+      sx3ZFrac = TMath::QuietNaN();
+      
+      sx3X = TMath::QuietNaN();
+      sx3Y = TMath::QuietNaN();
+      sx3Z = TMath::QuietNaN();
+
+      // for( int i = 0; i < 12; i++){
+      //   sx3Index[i] = -1;
+      // }
+  
+      reTheta = TMath::QuietNaN();
+      rePhi = TMath::QuietNaN();  
+
+      reTheta1 = TMath::QuietNaN();
+      rePhi1 = TMath::QuietNaN();
+
+      z0 = TMath::QuietNaN();
+
+    }
+    // -----------------------------------------------------------------------
+
+
+  }
+
+  // update TEve
+  pts->Reset();
+  for(size_t i=0;i<vx.size(); ++i) pts->SetNextPoint(vx[i], vy[i], vz[i]);
+  gEve->Redraw3D();
+  gSystem->ProcessEvents();
+
+  // write to tree
+  if(tout){ tout->Fill(); fout->Flush(); }
+
+  if(fout) fout->Close();
+  return 0;
+}
+
+// optional main to keep standalone build working
+#ifndef __CLING__
+int main(int argc, char** argv){
+  TApplication app("app",&argc,argv);
+  // if you want to import geometry here when running standalone:
+  // TGeoManager::Import("yourGeom.root");
+  Run(500, "sim_out.root");
+  return 0;
+}
+#endif
diff --git a/Armory/run_script.C b/Armory/run_script.C
new file mode 100644
index 0000000..e798df6
--- /dev/null
+++ b/Armory/run_script.C
@@ -0,0 +1,4 @@
+.L ANASEN_model.C
+.L anasenMS_root.cpp+
+ANASEN_model();
+Run(10);
diff --git a/Armory/test_run.C b/Armory/test_run.C
new file mode 100644
index 0000000..120b12d
--- /dev/null
+++ b/Armory/test_run.C
@@ -0,0 +1,7 @@
+.L ANASEN_model.C
+.L anasenMS_root.cpp+
+
+void test_run(){
+    ANASEN_model();
+    Run(10);
+}
diff --git a/Armory/vis_helpers.h b/Armory/vis_helpers.h
new file mode 100644
index 0000000..b6c01c8
--- /dev/null
+++ b/Armory/vis_helpers.h
@@ -0,0 +1,34 @@
+// vis_helpers.h  (or paste into anasenMS.cpp)
+#include <TEvePointSet.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <vector>
+#include <mutex>
+
+static TEvePointSet* gVisPts = nullptr;
+static std::mutex gVisMutex;
+
+// Call from your ROOT session after creating TEve objects:
+void SetVisPointSet(TEvePointSet* pts){ gVisPts = pts; }
+
+// Call this from your sim loop to update visualization and optionally write data:
+void PushEventAndRecord(const std::vector<double>& x,
+                        const std::vector<double>& y,
+                        const std::vector<double>& z,
+                        TTree* outTree = nullptr)
+{
+  if(outTree){
+    outTree->SetBranchAddress("x",(void*)&x);
+    outTree->SetBranchAddress("y",(void*)&y);
+    outTree->SetBranchAddress("z",(void*)&z);
+    outTree->Fill();
+    outTree->GetCurrentFile()->Flush();
+  }
+
+  if(!gVisPts) return;
+  std::lock_guard<std::mutex> lk(gVisMutex);
+  gVisPts->Reset();
+  for(size_t i=0;i<x.size(); ++i) gVisPts->SetNextPoint(x[i], y[i], z[i]);
+  gEve->Redraw3D();
+  gSystem->ProcessEvents();
+}
diff --git a/Armory/vis_inproc b/Armory/vis_inproc
new file mode 100755
index 0000000..390267e
Binary files /dev/null and b/Armory/vis_inproc differ
diff --git a/Armory/vis_inproc.cpp b/Armory/vis_inproc.cpp
new file mode 100644
index 0000000..29a6865
--- /dev/null
+++ b/Armory/vis_inproc.cpp
@@ -0,0 +1,65 @@
+//In-ROOT sim data
+
+// vis_inproc.cpp
+#include <TApplication.h>
+#include <TGeoManager.h>
+#include <TEveManager.h>
+#include <TEvePointSet.h>
+#include <TSystem.h>
+#include <TRandom3.h>
+#include <thread>
+#include <chrono>
+#include <vector>
+#include <iostream>
+
+void runSimulationAndUpdate(TEvePointSet* pts){
+  TRandom3 rnd(0);
+  for(int ev=0; ev<10000; ++ev){
+    pts->Reset();
+    int n = 100;
+    for(int i=0;i<n;++i){
+      double x = rnd.Uniform(-50,50);
+      double y = rnd.Uniform(-50,50);
+      double z = rnd.Uniform(-50,50);
+      pts->SetNextPoint(x,y,z);
+    }
+    gEve->Redraw3D();
+    gSystem->ProcessEvents();
+    std::this_thread::sleep_for(std::chrono::milliseconds(100));
+  }
+}
+
+void StartVis(const char* geomfile=nullptr){
+  int argc=0; char** argv=nullptr;
+  TApplication app("app",&argc,argv);
+  if(geomfile) TGeoManager::Import(geomfile);
+  TEveManager::Create();
+  TEvePointSet *pts = new TEvePointSet("hits");
+  gEve->AddElement(pts);
+  // runSimulationAndUpdate(pts);  // or leave update API to caller
+  app.Run();
+}
+
+int main(int argc, char** argv){
+  TApplication app("app",&argc,argv);
+
+  if(argc>1) TGeoManager::Import(argv[1]);
+
+  TEveManager::Create();
+
+  TEvePointSet *pts = new TEvePointSet("hits");
+  pts->SetMarkerStyle(20);
+  pts->SetMarkerSize(1.2);
+  pts->SetMarkerColor(kRed);
+  gEve->AddElement(pts);
+
+  // Option A: run simulation in same thread but yield to event loop inside the sim
+  runSimulationAndUpdate(pts);
+
+  // Option B: run sim in a separate thread (only if sim avoids ROOT globals)
+  // std::thread simThread(runSimulationAndUpdate, pts);
+  // simThread.detach();
+
+  app.Run();
+  return 0;
+}