1) pc-calibration macros added, 2) mild bookkeeping changes to handle multiple datasets with minimal changes. more changes to delta-t etc still pending

Armory/Armory/ANASEN_model.C

@@ -0,0 +1,121 @@
+#include "TGeoManager.h"
+#include "TGeoVolume.h"
+#include "TGeoBBox.h"
+#include "TCanvas.h"
+#include "TPolyMarker3D.h"
+#include "TPolyLine3D.h"
+
+#include "TMath.h"
+
+void ANASEN_model(int anodeID1 = -1, int anodeID2 = -1, int cathodeID1 = -1, int cathodeID2 = -1 ) {
+
+  // Create ROOT manager and master volume
+  TGeoManager *geom = new TGeoManager("Detector", "ANASEN");
+
+  //--- define some materials
+  TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
+  TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);
+  //--- define some media
+  TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
+  TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
+
+  //--- make the top container volume 
+  Double_t worldx = 200.; //mm
+  Double_t worldy = 200.; //mm
+  Double_t worldz = 200.; //mm
+  TGeoVolume *worldBox = geom->MakeBox("ROOT", Vacuum, worldx, worldy, worldz);
+  geom->SetTopVolume(worldBox);
+
+  //--- making axis
+  TGeoVolume *axisX = geom->MakeTube("axisX", Al, 0, 0.1, 5.);
+  axisX->SetLineColor(1);
+  worldBox->AddNode(axisX, 1, new TGeoCombiTrans(5, 0, 0., new TGeoRotation("rotA", 90., 90., 0.)));
+
+  TGeoVolume *axisY = geom->MakeTube("axisY", Al, 0, 0.1, 5.);
+  axisY->SetLineColor(1);
+  worldBox->AddNode(axisY, 1, new TGeoCombiTrans(0, 5, 0., new TGeoRotation("rotB", 0., 90., 0.)));
+
+  TGeoVolume *axisZ = geom->MakeTube("axisZ", Al, 0, 0.1, 5.);
+  axisZ->SetLineColor(1);
+  worldBox->AddNode(axisZ, 1, new TGeoTranslation(0, 0,  5));
+
+  //--- making ANASEN
+  const int nWire = 24;
+  const int wireShift = 3;
+  const int zLen = 300; //mm
+  const int radiusA = 38;
+  const int radiusC = 43;
+
+  //.......... convert to wire center dimensions
+  double dAngle = wireShift * TMath::TwoPi() / nWire;
+  double radiusAnew = radiusA * TMath::Cos( dAngle / 2.);
+  double wireALength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusA * TMath::Sin(dAngle/2),2) );
+  double wireATheta = TMath::ATan2( 2* radiusA * TMath::Sin( dAngle / 2.), zLen);
+
+  // printf("    dAngle : %f\n", dAngle);
+  // printf(" newRadius : %f\n", radiusAnew);
+  // printf("wireLength : %f\n", wireALength);
+  // printf("wire Theta : %f\n", wireATheta);
+
+  TGeoVolume *pcA = geom->MakeTube("tub1", Al, 0, 0.01, wireALength/2);
+  pcA->SetLineColor(4);  
+
+  for( int i = 0; i < nWire; i++){
+    if( anodeID2 >= 0 &&  (i < anodeID1 || i > anodeID2) ) continue;   
+    worldBox->AddNode(pcA, i+1, new TGeoCombiTrans( radiusAnew * TMath::Cos( TMath::TwoPi() / nWire *i + dAngle / 2), 
+                                                    radiusAnew * TMath::Sin( TMath::TwoPi() / nWire *i + dAngle / 2), 
+                                                    0, 
+                                                    new TGeoRotation("rot1", 360/ nWire * (i + wireShift/2.), wireATheta * 180/ TMath::Pi(), 0.)));
+
+  }
+
+  double radiusCnew = radiusC * TMath::Cos( dAngle / 2.);
+  double wireCLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusC * TMath::Sin(dAngle/2),2) );
+  double wireCTheta = TMath::ATan2( 2* radiusC * TMath::Sin( dAngle / 2.), zLen);
+
+  TGeoVolume *pcC = geom->MakeTube("tub2", Al, 0, 0.01, wireCLength/2);
+  pcC->SetLineColor(6);
+  for( int i = 0; i < nWire; i++){
+    if( cathodeID2 >= 0 && (i < cathodeID1 || i > cathodeID2) ) continue;   
+    worldBox->AddNode(pcC, i+1, new TGeoCombiTrans( radiusCnew * TMath::Cos( TMath::TwoPi() / nWire *i - dAngle/2), 
+                                                    radiusCnew * TMath::Sin( TMath::TwoPi() / nWire *i - dAngle/2), 
+                                                    0, 
+                                                    new TGeoRotation("rot1", 360/ nWire * (i - wireShift/2.), -wireCTheta * 180/ TMath::Pi(), 0.)));
+  }
+
+  const int nSX3 = 12;
+  const int sx3Radius = 88;
+  const int sx3Width = 40;
+  const int sx3Length = 75;
+  const int sx3Gap = 5;
+
+  TGeoVolume * sx3 = geom->MakeBox("box", Al, 0.1, sx3Width/2, sx3Length/2);
+  sx3->SetLineColor(kGreen+3);
+  for( int i = 0; i < nSX3; i++){
+    worldBox->AddNode(sx3, 2*i+1., new TGeoCombiTrans( sx3Radius * TMath::Cos( TMath::TwoPi() / nSX3 * (i + 0.5)), 
+                                                       sx3Radius * TMath::Sin( TMath::TwoPi() / nSX3 * (i + 0.5)), 
+                                                     sx3Length/2+sx3Gap, 
+                                                     new TGeoRotation("rot1", 360/nSX3 * (i + 0.5), 0., 0.)));
+
+    worldBox->AddNode(sx3, 2*i+2., new TGeoCombiTrans( sx3Radius * TMath::Cos( TMath::TwoPi() / nSX3 * (i + 0.5)), 
+                                                       sx3Radius * TMath::Sin( TMath::TwoPi() / nSX3 * (i + 0.5)), 
+                                                       -sx3Length/2-sx3Gap, 
+                                                       new TGeoRotation("rot1", 360/nSX3 * (i + 0.5), 0., 0.)));
+  }
+
+  const int qqqR1 = 50;
+  const int qqqR2 = 100;
+  TGeoVolume *qqq = geom->MakeTubs("qqq", Al, qqqR1, qqqR2, 0.5, 5, 85);
+  qqq->SetLineColor(7);
+  for( int i = 0; i < 4; i++){
+    worldBox->AddNode(qqq, i+1, new TGeoCombiTrans( 0, 
+                                                    0, 
+                                                    100, 
+                                                    new TGeoRotation("rot1", 360/4 * (i), 0., 0.)));
+  }
+
+  geom->CloseGeometry();
+
+  geom->SetVisLevel(4);
+  worldBox->Draw("ogle");
+}

Armory/Armory/ClassAnasen.h

@@ -0,0 +1,314 @@
+#ifndef ClassAnasen_h
+#define ClassAnasen_h
+
+#include <cstdio>
+#include <TMath.h>
+#include <TVector3.h>
+
+#include "TGeoManager.h"
+#include "TGeoVolume.h"
+#include "TGeoBBox.h"
+#include "TCanvas.h"
+#include "TPolyMarker3D.h"
+#include "TPolyLine3D.h"
+#include "TRandom.h"
+
+#include "ClassSX3.h"
+#include "ClassPW.h"
+
+class ANASEN{
+public:
+  ANASEN();
+  ~ANASEN();
+
+  void SetUncertainties(double sx3W, double sx3L, double anode, double cathode){
+    sigmaA = anode;
+    sigmaC = cathode;
+    sigmaW = sx3W;
+    sigmaL = sx3L;
+  }
+
+  void DrawTrack(TVector3 pos, TVector3 direction, bool drawEstimatedTrack = false);
+  void DrawDeducedTrack(TVector3 sx3Pos, int anodeID, int cathodeID);
+  void DrawAnasen(int anodeID1 = -1, 
+                  int anodeID2 = -1, 
+                  int cathodeID1 = -1, 
+                  int cathodeID2 = -1, 
+                  int sx3ID = -1, 
+                  bool DrawQQQ = false );
+
+
+
+  PW * GetPW() {return pw;}
+  SX3 * GetSX3() {return sx3;}
+
+private:
+
+  PW * pw;
+  SX3 * sx3;
+
+  double sigmaA, sigmaC; // pw
+  double sigmaW, sigmaL; // sx3
+
+  const float qqqR1 = 50;
+  const float qqqR2 = 100;
+  const float qqqZPos = 23 + 75 + 30; 
+
+  void CalGeometry();
+
+  TGeoManager *geom;
+  TGeoVolume *worldBox;
+
+  void Construct3DModel(int anodeID1 = -1, 
+                        int anodeID2 = -1, 
+                        int cathodeID1 = -1, 
+                        int cathodeID2 = -1, 
+                        int sx3ID = -1, 
+                        bool DrawQQQ = true);
+
+};
+
+//!==============================================
+inline ANASEN::ANASEN(){
+
+  pw = new PW();
+  sx3 = new SX3();
+
+  CalGeometry();
+
+  geom = nullptr;
+  worldBox = nullptr;
+
+}
+
+inline ANASEN::~ANASEN(){
+
+  delete geom;
+
+  delete pw;
+  delete sx3;
+
+}
+
+//!==============================================
+inline void ANASEN::CalGeometry(){
+
+  sx3->ConstructGeo();
+  pw->ConstructGeo();
+
+}
+
+inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1, int cathodeID2, int sx3ID, bool DrawQQQ ){
+
+  if( geom ) delete geom;
+
+  // Create ROOT manager and master volume
+  geom = new TGeoManager("Detector", "ANASEN");
+
+  //--- define some materials
+  TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
+  TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);
+  //--- define some media
+  TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
+  TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
+
+  //--- make the top container volume 
+  Double_t worldx = 200.; //mm
+  Double_t worldy = 200.; //mm
+  Double_t worldz = 200.; //mm
+  worldBox = geom->MakeBox("ROOT", Vacuum, worldx, worldy, worldz);
+  geom->SetTopVolume(worldBox);
+
+  //--- making axis
+  TGeoVolume *axisX = geom->MakeTube("axisX", Al, 0, 0.1, 5.);
+  axisX->SetLineColor(1);
+  worldBox->AddNode(axisX, 1, new TGeoCombiTrans(5, 0, 0., new TGeoRotation("rotA", 90., 90., 0.)));
+  TGeoVolume *axisY = geom->MakeTube("axisY", Al, 0, 0.1, 5.);
+  axisY->SetLineColor(1);
+  worldBox->AddNode(axisY, 1, new TGeoCombiTrans(0, 5, 0., new TGeoRotation("rotB", 0., 90., 0.)));
+  TGeoVolume *axisZ = geom->MakeTube("axisZ", Al, 0, 0.1, 5.);
+  axisZ->SetLineColor(1);
+  worldBox->AddNode(axisZ, 1, new TGeoTranslation(0, 0,  5));
+
+  //.......... convert to wire center dimensions
+  TGeoVolume *pcA = geom->MakeTube("tub1", Al, 0, 0.01, pw->GetAnodeLength()/2);
+  pcA->SetLineColor(4);  
+
+  int startID = 0;
+  int endID = pw->GetNumWire() - 1;
+
+  if( anodeID1 >= 0 && anodeID2 >= 0 ){
+    startID = anodeID1;
+    endID = anodeID2;
+    if( anodeID1 > anodeID2 ) {
+      endID = pw->GetNumWire() + anodeID2;
+    }
+  }
+
+  for( int i = startID; i <= endID; i++){
+    TVector3 a = pw->GetAnodneMid(i);
+    double wireTheta = pw->GetAnodeTheta(i) * TMath::RadToDeg();
+    double wirePhi = pw->GetAnodePhi(i) * TMath::RadToDeg() + 90;
+
+    worldBox->AddNode(pcA, i+1, new TGeoCombiTrans( a.X(), 
+                                                    a.Y(), 
+                                                    a.Z(), 
+                                                    new TGeoRotation("rot1", wirePhi, wireTheta, 0.)));
+  }
+
+  TGeoVolume *pcC = geom->MakeTube("tub2", Al, 0, 0.01, pw->GetCathodeLength()/2);
+  pcC->SetLineColor(6);
+
+  startID = 0;
+  endID = pw->GetNumWire() - 1;
+
+  if( cathodeID1 >= 0 && cathodeID2 >= 0 ){
+    startID = cathodeID1;
+    endID = cathodeID2;
+    if( cathodeID1 > cathodeID2 ) {
+      endID = pw->GetNumWire() + cathodeID2;
+    }
+  }
+
+  for( int i = startID; i <= endID; i++){
+    TVector3 a = pw->GetCathodneMid(i);
+    double wireTheta = pw->GetCathodeTheta(i) * TMath::RadToDeg();
+    double wirePhi = pw->GetCathodePhi(i) * TMath::RadToDeg() + 90;
+
+    worldBox->AddNode(pcC, i+1, new TGeoCombiTrans( a.X(), 
+                                                    a.Y(), 
+                                                    a.Z(), 
+                                                    new TGeoRotation("rot1", wirePhi , wireTheta, 0.)));
+  }
+
+  TGeoVolume * sx3Det = geom->MakeBox("box", Al, 0.1, sx3->GetWidth()/2, sx3->GetLength()/2);
+  sx3Det->SetLineColor(kGreen+3);
+
+  for( int i = 0; i < sx3->GetNumDet(); i++){
+    if( sx3ID != -1 && i != sx3ID ) continue;
+    TVector3 aUp = sx3->GetUpMid(i); // center of the SX3 upstream
+    TVector3 aDn = sx3->GetDnMid(i); // center of the SX3 Downstream
+    double phi = sx3->GetDetPhi(i) * TMath::RadToDeg() + 90;
+
+    worldBox->AddNode(sx3Det, 2*i+1., new TGeoCombiTrans( aUp.X(), 
+                                                          aUp.Y(), 
+                                                          aUp.Z(), 
+                                                          new TGeoRotation("rot1", phi, 0., 0.)));
+    worldBox->AddNode(sx3Det, 2*i+1., new TGeoCombiTrans( aDn.X(), 
+                                                          aDn.Y(), 
+                                                          aDn.Z(), 
+                                                          new TGeoRotation("rot1", phi, 0., 0.)));
+  }
+
+  if( DrawQQQ ){
+    TGeoVolume *qqq = geom->MakeTubs("qqq", Al, qqqR1, qqqR2, 0.5, 5, 85);
+    qqq->SetLineColor(7);
+    for( int i = 0; i < 4; i++){
+      worldBox->AddNode(qqq, i+1, new TGeoCombiTrans( 0, 
+                                                      0, 
+                                                      qqqZPos, 
+                                                      new TGeoRotation("rot1", 360/4 * (i), 0., 0.)));
+    }
+  }
+
+}
+
+//!============================================== Drawing functions
+inline void ANASEN::DrawAnasen(int anodeID1, int anodeID2, int cathodeID1, int cathodeID2, int sx3ID, bool DrawQQQ ){
+
+  Construct3DModel(anodeID1, anodeID2, cathodeID1, cathodeID2, sx3ID, DrawQQQ);
+
+  geom->CloseGeometry();
+  geom->SetVisLevel(4);
+  worldBox->Draw("ogle");
+
+}
+
+inline  void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstimatedTrack){
+
+  pw->FindWireID(pos, direction);
+  sx3->FindSX3Pos(pos, direction);
+
+  std::pair<short, short> wireID = pw->GetNearestID();
+  
+  Construct3DModel(wireID.first, wireID.first, wireID.second, wireID.second, -1, false);
+
+  double theta = direction.Theta() * TMath::RadToDeg();
+  double phi = direction.Phi()  * TMath::RadToDeg();
+  // printf("Theta, Phi = %.2f %.2f \n", theta, phi);
+  // pos.Print();
+  TGeoVolume * Track = geom->MakeTube("track", 0, 0, 0.1, 150.);
+  Track->SetLineColor(kRed);
+  worldBox->AddNode(Track, 1, new TGeoCombiTrans( pos.X(), pos.Y(), pos.Z(), new TGeoRotation("rotA", phi +  90, theta, 0.)));
+
+  TGeoVolume * startPos = geom->MakeSphere("startPos", 0, 0, 3);
+  startPos->SetLineColor(kBlack);
+  worldBox->AddNode(startPos, 3, new TGeoCombiTrans( pos.X(), pos.Y(), pos.Z(), new TGeoRotation("rotA", 0, 0, 0.)));
+
+  if( sx3->GetID() >= 0 ){
+    //TVector3 hitPos = sx3->GetHitPos();
+    TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaW, sigmaL);
+
+    TGeoVolume * hit = geom->MakeSphere("hitpos", 0, 0, 3);
+    hit->SetLineColor(kRed);
+    worldBox->AddNode(hit, 2, new TGeoCombiTrans( hitPos.X(), hitPos.Y(), hitPos.Z(), new TGeoRotation("rotA", 0, 0, 0.)));
+
+    if( drawEstimatedTrack ){
+
+      {//===== simple
+        pw->CalTrack(hitPos, wireID.first, wireID.second, true);
+
+        double thetaDeduce = pw->GetTrackTheta() * TMath::RadToDeg();
+        double phiDeduce = pw->GetTrackPhi()  * TMath::RadToDeg();
+
+        TGeoVolume * trackDeduce = geom->MakeTube("trackDeduce", 0, 0, 0.1, 100.);
+        trackDeduce->SetLineColor(kOrange);
+        worldBox->AddNode(trackDeduce, 1, new TGeoCombiTrans( hitPos.X(), hitPos.Y(), hitPos.Z(), new TGeoRotation("rotA", phiDeduce +  90, thetaDeduce, 0.)));
+      }
+
+      {//===== complicated 
+        PWHitInfo hitInfo = pw->GetHitInfo();
+        pw->CalTrack2(hitPos, hitInfo, sigmaA, sigmaC, true);
+
+        double thetaDeduce = pw->GetTrackTheta() * TMath::RadToDeg();
+        double phiDeduce = pw->GetTrackPhi()  * TMath::RadToDeg();
+
+        TGeoVolume * trackDeduce2 = geom->MakeTube("trackDeduce2", 0, 0, 0.1, 100.);
+        trackDeduce2->SetLineColor(kGreen);
+        worldBox->AddNode(trackDeduce2, 1, new TGeoCombiTrans( hitPos.X(), hitPos.Y(), hitPos.Z(), new TGeoRotation("rotA", phiDeduce +  90, thetaDeduce, 0.)));
+      }
+    }
+
+  }
+  
+  geom->CloseGeometry();
+  geom->SetVisLevel(4);
+  worldBox->Draw("ogle");
+
+}
+
+inline  void ANASEN::DrawDeducedTrack(TVector3 sx3Pos, int anodeID, int cathodeID){
+
+  pw->CalTrack(sx3Pos, anodeID, cathodeID);
+
+  Construct3DModel(anodeID, anodeID, cathodeID, cathodeID, -1, false);
+
+  double theta = pw->GetTrackTheta() * TMath::RadToDeg();
+  double phi = pw->GetTrackPhi()  * TMath::RadToDeg();
+
+  TGeoVolume * Track = geom->MakeTube("axisX", 0, 0, 0.1, 100.);
+  Track->SetLineColor(kRed);
+  worldBox->AddNode(Track, 1, new TGeoCombiTrans( sx3Pos.X(), sx3Pos.Y(), sx3Pos.Z(), new TGeoRotation("rotA", phi +  90, theta, 0.)));
+
+  TGeoVolume * hit = geom->MakeSphere("hitpos", 0, 0, 3);
+  hit->SetLineColor(kRed);
+  worldBox->AddNode(hit, 2, new TGeoCombiTrans( sx3Pos.X(), sx3Pos.Y(), sx3Pos.Z(), new TGeoRotation("rotA", 0, 0, 0.)));
+
+  geom->CloseGeometry();
+  geom->SetVisLevel(4);
+  worldBox->Draw("ogle");
+
+}
+
+
+#endif 

Armory/Armory/ClassDet.h

@@ -0,0 +1,64 @@
+#ifndef ClassDet_h
+#define ClassDet_h
+
+#include <cstdio>
+
+#define MAXMULTI 1000
+
+class Det{
+public:
+  Det(): multi(0) {Clear(); }
+
+  unsigned short       multi; // max 65535
+  unsigned short       id[MAXMULTI];
+  unsigned short       ch[MAXMULTI];
+  unsigned short       e[MAXMULTI];
+  unsigned long long   t[MAXMULTI];
+
+  unsigned short         sn[MAXMULTI];
+  unsigned short     digiCh[MAXMULTI];
+
+  unsigned short  index[MAXMULTI]; // id * nCh + ch;
+  bool  used[MAXMULTI];
+
+  void Clear(){
+    multi = 0;
+    for( int i = 0; i < MAXMULTI; i++){
+      id[i] = 0;
+      ch[i] = 0;
+       e[i] = 0;
+       t[i] = 0;
+       index[i] = 0;
+       sn[i] = 0;
+       digiCh[i] = 0;
+       used[i] = false;
+    }
+  }
+
+  void Print(){
+    printf("=============================== multi : %u\n", multi);
+    for( int i = 0; i < multi; i++) {
+      printf(" %3d | %2d-%-2d(%5d) %5u %15llu \n", i, id[i], ch[i], index[i], e[i], t[i]);
+    }
+  }
+
+  void SetDetDimension(unsigned short maxID, unsigned maxCh){
+    nID = maxID;
+    nCh = maxCh;
+  }
+
+  void CalIndex(){
+    for( int i = 0; i < multi; i++){
+      index[i] = id[i] * nCh + ch[i] ;
+    }
+  }
+
+private:
+
+  unsigned short nID;
+  unsigned short nCh;
+
+};
+
+
+#endif 

Armory/Armory/ClassPC1An.h

@@ -0,0 +1,246 @@
+#ifndef ClassPC_h
+#define ClassPC_h
+
+#include <cstdio>
+#include <TMath.h>
+#include <TVector3.h>
+#include <TRandom.h>
+
+struct PCHit_1An{
+  std::pair<short, short>   nearestWire; // anode, cathode
+  std::pair<double, double> nearestDist; // anode, cathode
+
+  short   nextNearestWire; // cathode
+  double nextNearestDist; // cathode
+
+  void Clear(){
+    nearestWire.first = -1;
+    nearestWire.second = -1;
+    nearestDist.first = 999999999;
+    nearestDist.second = 999999999;
+    nextNearestWire= -1;
+    nextNearestDist = 999999999;
+  }
+};
+
+//!########################################################
+class PC{ // proportional wire
+public: 
+  PC(){ ClearHitInfo();};
+  ~PC(){};
+  
+  PCHit_1An GetHitInfo() const {return hitInfo;}
+  std::pair<short, short>   GetNearestID()          const {return hitInfo.nearestWire;}
+  std::pair<double, double> GetNearestDistance()    const {return hitInfo.nearestDist;}
+  short   Get2ndNearestID()       const {return hitInfo.nextNearestWire;}
+  double Get2ndNearestDistance() const {return hitInfo.nextNearestDist;}
+
+  TVector3 GetTrackPos() const {return trackPos;}
+  TVector3 GetTrackVec() const {return trackVec;}
+  double GetTrackTheta() const {return trackVec.Theta();}
+  double GetTrackPhi()   const {return trackVec.Phi();}
+  double GetZ0();
+
+  int GetNumWire() const {return nWire;}
+  double GetDeltaAngle() const {return dAngle;}
+  double GetAnodeLength() const {return anodeLength;}
+  double GetCathodeLength() const {return cathodeLength;}
+  TVector3 GetAnodeDn(short id) const {return An[id].first;}
+  TVector3 GetAnodeUp(short id) const {return An[id].second;}
+  TVector3 GetCathodeDn(short id) const {return Ca[id].first;}
+  TVector3 GetCathodeUp(short id) const {return Ca[id].second;}
+
+  TVector3 GetAnodneMid(short id)  const {return (An[id].first + An[id].second) * 0.5; }
+  double   GetAnodeTheta(short id) const {return (An[id].first - An[id].second).Theta();}
+  double   GetAnodePhi(short id)   const {return (An[id].first - An[id].second).Phi();}
+
+  TVector3 GetCathodneMid(short id)  const {return (Ca[id].first + Ca[id].second) * 0.5; }
+  double   GetCathodeTheta(short id) const {return (Ca[id].first - Ca[id].second).Theta();}
+  double   GetCathodePhi(short id)   const {return (Ca[id].first - Ca[id].second).Phi();}
+
+  void ClearHitInfo();
+  void ConstructGeo();
+  void FindWireID(TVector3 pos, TVector3 direction, bool verbose = false);
+  void CalTrack3(TVector3 sx3Pos, PCHit_1An hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false);
+
+  void Print(){
+    printf("     The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first, 
+                                                                         hitInfo.nearestDist.first, 
+                                                                         hitInfo.nearestWire.second, 
+                                                                         hitInfo.nearestDist.second);
+
+    printf(" The 2nd nearest Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire, 
+                                                     hitInfo.nextNearestDist);
+  }
+
+private:
+
+  // PCHitInfo hitInfo;
+  PCHit_1An hitInfo;
+
+  TVector3 trackPos;
+  TVector3 trackVec;
+
+  const int nWire = 24;
+  const int wireShift = 3;
+  const float zLen = 380; //mm
+  const float radiusA = 37;
+  const float radiusC = 43;
+
+  double dAngle;
+  double anodeLength;
+  double cathodeLength;
+
+  std::vector<std::pair<TVector3,TVector3>> An; // the anode wire position vector in space 
+  std::vector<std::pair<TVector3,TVector3>> Ca; // the cathode wire position vector in space 
+
+  double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2){
+    TVector3 na = a1 - a2;
+    TVector3 nb = b1 - b2;
+    TVector3 nd = (na.Cross(nb)).Unit();
+    return TMath::Abs(nd.Dot(a1-b2));
+  } 
+
+};
+
+inline void PC::ClearHitInfo(){
+  hitInfo.Clear();
+}
+
+inline void PC::ConstructGeo(){
+
+  An.clear();
+  Ca.clear();
+
+  std::pair<TVector3, TVector3> p1; // anode
+  std::pair<TVector3, TVector3> q1; // cathode
+
+  //anode and cathode start at pos-Y axis and count in right-Hand
+  //anode wire shift is right-hand.
+  //cathode wire shift is left-hand.
+
+  for(int i = 0; i < nWire; i++ ){
+    // Anode rotate right-hand
+    p1.first.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i)  + TMath::PiOver2()),
+                    radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i)  + TMath::PiOver2()),
+                    zLen/2);
+    p1.second.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
+                      radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
+                    -zLen/2);
+    An.push_back(p1);
+
+    // Cathod rotate left-hand
+    q1.first.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i)  + TMath::PiOver2()),
+                    radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i)  + TMath::PiOver2()),
+                    zLen/2);
+    q1.second.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
+                      radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
+                      -zLen/2);
+    Ca.push_back(q1);
+  }
+
+  dAngle = wireShift * TMath::TwoPi() / nWire;
+  anodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusA * TMath::Sin(dAngle/2),2) );
+  cathodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusC * TMath::Sin(dAngle/2),2) );
+}
+
+inline void PC::FindWireID(TVector3 pos, TVector3 direction, bool verbose ){
+
+  hitInfo.Clear();
+  double phi = direction.Phi();
+
+  for( int i = 0; i < nWire; i++){
+
+    double disA = 99999999;
+    double phiS = An[i].first.Phi()  - TMath::PiOver4();
+    double phiL = An[i].second.Phi() + TMath::PiOver4();
+    // printf("A%2d: %f %f | %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg(), phi * TMath::RadToDeg());
+    if( phi > 0 && phiS > phiL )  phiL = phiL + TMath::TwoPi();
+    if( phi < 0 && phiS > phiL )  phiS = phiS - TMath::TwoPi();
+
+    if( phiS < phi && phi < phiL) {
+      disA = Distance( pos, pos + direction, An[i].first, An[i].second);
+      if( disA < hitInfo.nearestDist.first ){
+        hitInfo.nearestDist.first = disA;
+        hitInfo.nearestWire.first = i;
+      }
+    }
+
+    double disC = 99999999;
+    phiS = Ca[i].second.Phi()- TMath::PiOver4();
+    phiL = Ca[i].first.Phi() + TMath::PiOver4();
+    // printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
+    if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi();
+    if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi();
+
+    if(phiS < phi && phi < phiL) {
+      disC = Distance( pos, pos + direction, Ca[i].first, Ca[i].second);
+      if( disC < hitInfo.nearestDist.second ){
+        hitInfo.nearestDist.second = disC;
+        hitInfo.nearestWire.second = i;
+      }
+    }
+
+    if(verbose) printf(" %2d | %8.2f, %8.2f\n", i, disA, disC);
+  }
+
+  short cathode1 = hitInfo.nearestWire.second;
+  short ccc1 = cathode1 - 1; if( ccc1 < 0 ) ccc1 += nWire;
+  short ccc2 = (cathode1 + 1) % nWire; 
+
+   double haha1 = Distance( pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second);
+   double haha2 = Distance( pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second);
+  if( haha1 < haha2){
+    hitInfo.nextNearestWire = ccc1;
+    hitInfo.nextNearestDist = haha1;
+  }else{
+    hitInfo.nextNearestWire = ccc2;
+    hitInfo.nextNearestDist= haha2;
+  }
+
+  if( verbose ) Print();
+}
+
+
+inline void PC::CalTrack3(TVector3 sx3Pos, PCHit_1An hitInfo, double sigmaA, double sigmaC, bool verbose){
+
+  trackPos = sx3Pos;
+
+  double p1 = TMath::Abs(hitInfo.nearestDist.first + gRandom->Gaus(0, sigmaA));
+  short anodeID1 = hitInfo.nearestWire.first;
+
+  double q1 = TMath::Abs(hitInfo.nearestDist.second + gRandom->Gaus(0, sigmaC));
+  double q2 = TMath::Abs(hitInfo.nextNearestDist+ gRandom->Gaus(0, sigmaC));
+  double fracC = q1 / (q1 + q2);
+  short cathodeID1 = hitInfo.nearestWire.second;
+  short cathodeID2 = hitInfo.nextNearestWire;
+  TVector3 shiftC1 = (Ca[cathodeID2].first - Ca[cathodeID1].first) * fracC;
+  TVector3 shiftC2 = (Ca[cathodeID2].second - Ca[cathodeID1].second) * fracC;
+
+  TVector3 a1 = An[anodeID1].first;
+
+  TVector3 c1 = Ca[cathodeID1].first + shiftC1;
+  TVector3 c2 = Ca[cathodeID1].second + shiftC2;
+
+  TVector3 n1 = (sx3Pos - a1).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("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg()); 
+
+}
+
+inline double PC::GetZ0(){
+
+  double x = trackPos.X();
+  double y = trackPos.Y();
+  double rho = TMath::Sqrt(x*x + y*y);
+  double theta = trackVec.Theta();
+  
+  return trackPos.Z() - rho / TMath::Tan(theta);
+
+}
+
+#endif 

Armory/Armory/ClassPW.h

@@ -0,0 +1,315 @@
+#ifndef ClassPW_h
+#define ClassPW_h
+
+#include <cstdio>
+#include <TMath.h>
+#include <TVector3.h>
+#include <TRandom.h>
+
+struct PWHitInfo
+{
+  std::pair<short, short> nearestWire;   // anode, cathode
+  std::pair<double, double> nearestDist; // anode, cathode
+
+  std::pair<short, short> nextNearestWire;   // anode, cathode
+  std::pair<double, double> nextNearestDist; // anode, cathode
+
+  void Clear()
+  {
+    nearestWire.first = -1;
+    nearestWire.second = -1;
+    nearestDist.first = 999999999;
+    nearestDist.second = 999999999;
+    nextNearestWire.first = -1;
+    nextNearestWire.second = -1;
+    nextNearestDist.first = 999999999;
+    nextNearestDist.second = 999999999;
+  }
+};
+
+struct Coord
+{
+  float x, y, z;
+  Coord() : x(0), y(0), z(0) {}
+  Coord(const TVector3 &vec)
+  {
+    x = vec.X(); // TVector3's X() returns the x-coordinate
+    y = vec.Y(); // TVector3's Y() returns the y-coordinate
+    z = vec.Z(); // TVector3's Z() returns the z-coordinate
+  }
+};
+
+//! ########################################################
+class PW
+{ // proportional wire
+public:
+  PW() { ClearHitInfo(); };
+  ~PW() {};
+
+  PWHitInfo GetHitInfo() const { return hitInfo; }
+  std::pair<short, short> GetNearestID() const { return hitInfo.nearestWire; }
+  std::pair<double, double> GetNearestDistance() const { return hitInfo.nearestDist; }
+  std::pair<short, short> Get2ndNearestID() const { return hitInfo.nextNearestWire; }
+  std::pair<double, double> Get2ndNearestDistance() const { return hitInfo.nextNearestDist; }
+
+  std::vector<std::pair<TVector3, TVector3>> An; // the anode wire position vector in space
+  std::vector<std::pair<TVector3, TVector3>> Ca; // the cathode wire position vector in space
+
+  TVector3 GetTrackPos() const { return trackPos; }
+  TVector3 GetTrackVec() const { return trackVec; }
+  double GetTrackTheta() const { return trackVec.Theta(); }
+  double GetTrackPhi() const { return trackVec.Phi(); }
+  double GetZ0();
+
+  int GetNumWire() const { return nWire; }
+  double GetDeltaAngle() const { return dAngle; }
+  double GetAnodeLength() const { return anodeLength; }
+  double GetCathodeLength() const { return cathodeLength; }
+  TVector3 GetAnodeDn(short id) const { return An[id].first; }
+  TVector3 GetAnodeUp(short id) const { return An[id].second; }
+  TVector3 GetCathodeDn(short id) const { return Ca[id].first; }
+  TVector3 GetCathodeUp(short id) const { return Ca[id].second; }
+
+  TVector3 GetAnodneMid(short id) const { return (An[id].first + An[id].second) * 0.5; }
+  double GetAnodeTheta(short id) const { return (An[id].first - An[id].second).Theta(); }
+  double GetAnodePhi(short id) const { return (An[id].first - An[id].second).Phi(); }
+
+  TVector3 GetCathodneMid(short id) const { return (Ca[id].first + Ca[id].second) * 0.5; }
+  double GetCathodeTheta(short id) const { return (Ca[id].first - Ca[id].second).Theta(); }
+  double GetCathodePhi(short id) const { return (Ca[id].first - Ca[id].second).Phi(); }
+
+  void ClearHitInfo();
+  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 Print()
+  {
+    printf("     The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first,
+           hitInfo.nearestDist.first,
+           hitInfo.nearestWire.second,
+           hitInfo.nearestDist.second);
+
+    printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire.first,
+           hitInfo.nextNearestDist.first,
+           hitInfo.nextNearestWire.second,
+           hitInfo.nextNearestDist.second);
+  }
+
+private:
+  PWHitInfo hitInfo;
+
+  TVector3 trackPos;
+  TVector3 trackVec;
+
+  const int nWire = 24;
+  const int wireShift = 3;
+  const float zLen = 380; // mm
+  const float radiusA = 37;
+  const float radiusC = 43;
+
+  double dAngle;
+  double anodeLength;
+  double cathodeLength;
+
+  // std::vector<std::pair<TVector3, TVector3>> An; // the anode wire position vector in space
+  // std::vector<std::pair<TVector3, TVector3>> Ca; // the cathode wire position vector in space
+
+  double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2)
+  {
+    TVector3 na = a1 - a2;
+    TVector3 nb = b1 - b2;
+    TVector3 nd = (na.Cross(nb)).Unit();
+    return TMath::Abs(nd.Dot(a1 - b2));
+  }
+};
+
+inline void PW::ClearHitInfo()
+{
+  hitInfo.Clear();
+}
+
+inline void PW::ConstructGeo()
+{
+
+  An.clear();
+  Ca.clear();
+
+  std::pair<TVector3, TVector3> p1; // anode
+  std::pair<TVector3, TVector3> q1; // cathode
+
+  // anode and cathode start at pos-Y axis and count in right-Hand
+  // anode wire shift is right-hand.
+  // cathode wire shift is left-hand.
+
+  for (int i = 0; i < nWire; i++)
+  {
+    // Anode rotate right-hand
+    p1.first.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
+                    radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
+                    zLen / 2);
+    p1.second.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
+                     radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
+                     -zLen / 2);
+    An.push_back(p1);
+
+    // Cathod rotate left-hand with the 3 wire offset accounted for (+1 from the calculated offset from the PC coincidence spectrum)
+    q1.first.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()),
+                    radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()),
+                    zLen / 2);
+    q1.second.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i + 1) + TMath::PiOver2()),
+                     radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i + 1) + TMath::PiOver2()),
+                     -zLen / 2);
+    Ca.push_back(q1);
+  }
+  // correcting for the fact that the order of the cathode wires is reversed
+  std::reverse(Ca.begin(), Ca.end());
+  // adjusting for the 3 wire offset, the rbegin and rend are used as the rotation of the wires is done in the opposite direction i.e. 1,2,3 -> 3,1,2
+  // NOT NECESSARY ANY MORE, HAS BEEN IMCORPORATED INTO THE WIREOFFSET IN THE BEGINNING
+  // std::rotate(Ca.rbegin(), Ca.rbegin() + 4, Ca.rend());
+
+  dAngle = wireShift * TMath::TwoPi() / nWire;
+  anodeLength = TMath::Sqrt(zLen * zLen + TMath::Power(2 * radiusA * TMath::Sin(dAngle / 2), 2));
+  cathodeLength = TMath::Sqrt(zLen * zLen + TMath::Power(2 * radiusC * TMath::Sin(dAngle / 2), 2));
+}
+
+inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
+{
+
+  hitInfo.Clear();
+  double phi = direction.Phi();
+
+  for (int i = 0; i < nWire; i++)
+  {
+
+    double disA = 99999999;
+    double phiS = An[i].first.Phi() - TMath::PiOver4();
+    double phiL = An[i].second.Phi() + TMath::PiOver4();
+    // printf("A%2d: %f %f | %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg(), phi * TMath::RadToDeg());
+    if (phi > 0 && phiS > phiL)
+      phiL = phiL + TMath::TwoPi();
+    if (phi < 0 && phiS > phiL)
+      phiS = phiS - TMath::TwoPi();
+
+    if (phiS < phi && phi < phiL)
+    {
+      disA = Distance(pos, pos + direction, An[i].first, An[i].second);
+      if (disA < hitInfo.nearestDist.first)
+      {
+        hitInfo.nearestDist.first = disA;
+        hitInfo.nearestWire.first = i;
+      }
+    }
+
+    double disC = 99999999;
+    phiS = Ca[i].second.Phi() - TMath::PiOver4();
+    phiL = Ca[i].first.Phi() + TMath::PiOver4();
+    // printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
+    if (phi > 0 && phiS > phiL)
+      phiL = phiL + TMath::TwoPi();
+    if (phi < 0 && phiS > phiL)
+      phiS = phiS - TMath::TwoPi();
+
+    if (phiS < phi && phi < phiL)
+    {
+      disC = Distance(pos, pos + direction, Ca[i].first, Ca[i].second);
+      if (disC < hitInfo.nearestDist.second)
+      {
+        hitInfo.nearestDist.second = disC;
+        hitInfo.nearestWire.second = i;
+      }
+    }
+
+    if (verbose)
+      printf(" %2d | %8.2f, %8.2f\n", i, disA, disC);
+  }
+
+  //==== find the 2nd nearest wire
+  short anode1 = hitInfo.nearestWire.first;
+  short aaa1 = anode1 - 1;
+  if (aaa1 < 0)
+    aaa1 += nWire;
+  short aaa2 = (anode1 + 1) % nWire;
+
+  double haha1 = Distance(pos, pos + direction, An[aaa1].first, An[aaa1].second);
+  double haha2 = Distance(pos, pos + direction, An[aaa2].first, An[aaa2].second);
+  if (haha1 < haha2)
+  {
+    hitInfo.nextNearestWire.first = aaa1;
+    hitInfo.nextNearestDist.first = haha1;
+  }
+  else
+  {
+    hitInfo.nextNearestWire.first = aaa2;
+    hitInfo.nextNearestDist.first = haha2;
+  }
+
+  short cathode1 = hitInfo.nearestWire.second;
+  short ccc1 = cathode1 - 1;
+  if (ccc1 < 0)
+    ccc1 += nWire;
+  short ccc2 = (cathode1 + 1) % nWire;
+
+  haha1 = Distance(pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second);
+  haha2 = Distance(pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second);
+  if (haha1 < haha2)
+  {
+    hitInfo.nextNearestWire.second = ccc1;
+    hitInfo.nextNearestDist.second = haha1;
+  }
+  else
+  {
+    hitInfo.nextNearestWire.second = ccc2;
+    hitInfo.nextNearestDist.second = haha2;
+  }
+
+  if (verbose)
+    Print();
+}
+
+inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose)
+{
+
+  trackPos = sx3Pos;
+
+  TVector3 n1 = (An[anodeID].first - An[anodeID].second).Cross((sx3Pos - An[anodeID].second)).Unit();
+  TVector3 n2 = (Ca[cathodeID].first - Ca[cathodeID].second).Cross((sx3Pos - Ca[cathodeID].second)).Unit();
+
+  // if the handiness of anode and cathode revered, it should be n2 cross n1
+  trackVec = (n2.Cross(n1)).Unit();
+
+  if (verbose)
+    printf("Theta, Phi = %f, %f \n", trackVec.Theta() * TMath::RadToDeg(), trackVec.Phi() * TMath::RadToDeg());
+}
+
+
+inline void PW::CalTrack2(TVector3 siPos, TVector3 anodeInt, 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);
+  }
+
+  if (verbose)
+    printf("X slope = %f and Y slope = %f \n", mx, my);
+}
+
+inline double PW::GetZ0()
+{
+
+  double x = trackPos.X();
+  double y = trackPos.Y();
+  double rho = TMath::Sqrt(x * x + y * y);
+  double theta = trackVec.Theta();
+
+  return trackVec.Z();
+}
+
+#endif

Armory/Armory/ClassSX3.h

@@ -0,0 +1,257 @@
+#ifndef ClassSX3_h
+#define ClassSX3_h
+
+#include <cstdio>
+#include <TMath.h>
+#include <TVector3.h>
+#include <TRandom.h>
+
+class SX3{
+public:
+  SX3(){Clear();};
+  ~SX3(){}
+  
+  short GetID() const {return id;}
+  short GetChUp() const {return chUp;}
+  short GetChDn() const {return chDn;}
+  short GetChBk() const {return chBk;}
+
+  TVector3 GetHitPos() const {return hitPos;}
+  TVector3 GetHitPosWithSigma(double sigmaY_mm, double sigmaZ_mm);
+
+  double GetZFrac() const {return zFrac;} // range from -0.5 to 0.5
+
+  void Clear();
+  void ConstructGeo();
+  void FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose = false);
+  void CalSX3Pos(unsigned short ID, unsigned short chUp, unsigned short chDown, unsigned short chBack, float eUp, float eDown);
+
+  double GetNumDet() const {return numDet;}
+  double GetWidth()  const {return width;}
+  double GetLength() const {return length;}
+  TVector3 GetDnL(short id) const {return SDn[id].first; } // lower strip ID
+  TVector3 GetDnH(short id) const {return SDn[id].second; } // higher strip ID
+  TVector3 GetUpL(short id) const {return SUp[id].first; } // lower strip ID
+  TVector3 GetUpH(short id) const {return SUp[id].second; } // higher strip ID
+
+  TVector3 GetDnMid(short id) const { return (SDn[id].first + SDn[id].second)*0.5;}
+  TVector3 GetUpMid(short id) const { return (SUp[id].first + SUp[id].second)*0.5;}
+
+  double GetDetPhi(short id) const { return (SUp[id].second - SUp[id].first).Phi();}
+
+  void Print(){ 
+    if( id == -1 ){
+      printf("Did not hit any SX3.\n"); 
+    }else{
+      printf("ID: %d, U,D,B: %d %d %d| zFrac : %.2f\n", id, chUp, chDn, chBk, zFrac); 
+      printf("Hit Pos: %.2f, %.2f, %.2f\n", hitPos.X(), hitPos.Y(), hitPos.Z()); 
+    }
+  }
+
+  // void CalZFrac(){
+  //   zFrac = (eUp - eDn)/(eUp + eDn); 
+  // }
+
+private:
+
+  const int numDet = 12;
+  const float radius = 88;
+  const float width = 40;
+  const float length = 75;
+  const float gap = 46;
+
+  short id; // -1 when no hit
+  short chUp;
+  short chDn;
+  short chBk;
+
+  double zFrac; // from +1 (downstream) to -1 (upstream)
+
+  double eUp;
+  double eDn;
+  double eBk;
+
+  TVector3 hitPos;
+
+  std::vector<std::pair<TVector3,TVector3>> SDn; // coners of the SX3  0-11, z = mid point
+  std::vector<std::pair<TVector3,TVector3>> SUp; // coners of the SX3 12-23, z = mid point
+  std::vector<TVector3> SNorml; // normal of the SX3 (outward)
+
+
+  std::pair<double, double> Intersect(TVector3 p1, TVector3 p2, TVector3 q1, TVector3 q2, bool verbose){
+
+    //see https://nukephysik101.wordpress.com/2023/12/30/intersect-between-2-line-segments/
+    //zero all z-component
+    TVector3 a0 = p1; a0.SetZ(0);
+    TVector3 a1 = p2; a1.SetZ(0);
+
+    TVector3 b0 = q1; b0.SetZ(0);
+    TVector3 b1 = q2; b1.SetZ(0);
+
+    double A = ((b0-b1).Cross(a0-a1)).Mag();
+
+    double h = ((b0-a0).Cross(b1-a0)).Z()/ A;
+    double k = ((a1-b0).Cross(a0-b0)).Z()/ A;
+
+    if( verbose ) printf(" ----h, k : %f, %f\n", h, k);
+
+    return std::pair<double,double>(h,k);
+  } 
+
+};
+
+inline void SX3::Clear(){
+  id = -1;
+  chUp = -1;
+  chDn = -1;
+  chBk = -1;
+  zFrac = TMath::QuietNaN();
+
+  eUp = TMath::QuietNaN();
+  eDn = TMath::QuietNaN();
+  eBk = TMath::QuietNaN();
+
+  SDn.clear();
+  SUp.clear();
+}
+
+inline void SX3::ConstructGeo(){
+  TVector3 sa, sb, sc, sn;
+
+  for(int i = 0; i < numDet; i++){
+    sa.SetXYZ( radius, -width/2, gap/2 + length/2 );
+    sb.SetXYZ( radius,  width/2, gap/2 + length/2 );
+
+    double rot = TMath::TwoPi() / numDet * (i - 0.5) - TMath::PiOver2();
+
+    sa.RotateZ( rot );
+    sb.RotateZ( rot );
+    SDn.push_back(std::pair<TVector3,TVector3>(sa,sb));
+
+    sc.SetXYZ( radius, -width/2, gap/2 );
+    sc.RotateZ( rot );
+
+    sn = ((sc-sa).Cross(sb-sa)).Unit();
+    SNorml.push_back(sn);
+
+    sa.SetXYZ( radius, -width/2, -gap/2 - length/2 );
+    sb.SetXYZ( radius,  width/2, -gap/2 - length/2 );
+
+    sa.RotateZ( rot );
+    sb.RotateZ( rot );
+    SUp.push_back(std::pair<TVector3,TVector3>(sa,sb));
+  }
+}
+
+
+inline void SX3::FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose){
+
+  id = -1;
+  for( int i = 0 ; i < numDet; i++){
+
+    if(verbose) printf(" %d ", i);
+    std::pair<double, double> frac = Intersect( pos, pos + direction, SDn[i].first, SDn[i].second, verbose);
+
+
+    if( frac.second < 0 || frac.second > 1 ) continue;
+    hitPos = pos + frac.first * direction;
+
+    double dis = hitPos.Dot(SNorml[i]);
+
+    if(verbose) {
+      printf("reduced distance : %f\n", dis);
+      printf(" %d*", (i+1)%numDet); 
+      Intersect( pos, pos + direction, SDn[(i+1)%numDet].first, SDn[(i+1)%numDet].second, verbose);
+    }
+
+    if( TMath::Abs(dis - radius) > 0.1 ) continue;
+
+    chDn = 2 * TMath::Floor(frac.second * 4);
+    chUp = chDn + 1;
+
+    double zPos = hitPos.Z();
+    if( (gap/2 < zPos && zPos < gap/2 + length ) || (-gap/2 - length < zPos && zPos < -gap/2 )  ){
+
+      id = zPos > 0 ? i : i + 12;
+
+      zFrac = zPos > 0 ?  (zPos - gap/2. - length/2.)/length : (zPos - ( - gap/2. - length/2.) )/length ;
+
+      chBk = TMath::Floor( (zFrac + 0.5) * 4 ) + 8;
+
+      if( verbose) Print();
+      return ;
+
+    }else{
+      if( verbose ) printf(" zPos out of sensitive region\n");
+    }
+  }
+
+  if( verbose) Print();
+}
+
+inline TVector3 SX3::GetHitPosWithSigma(double sigmaY_mm, double sigmaZ_mm){
+
+  double phi = SNorml[id%numDet].Phi();
+
+  TVector3 haha = hitPos;
+  haha.RotateZ(-phi);
+
+  double y = haha.Y() + gRandom->Gaus(0, sigmaY_mm);
+  if( sigmaY_mm < 0 ){
+    double deltaW = width/4;
+    y = TMath::Floor((haha.Y()-deltaW)/deltaW)*deltaW + deltaW*1.5; // when ever land on each strip, set the position to be center of the strip.
+    if( y >= 25 ) y = 15;
+  }
+
+  double z = haha.Z() + gRandom->Gaus(0, sigmaZ_mm);
+  if( sigmaZ_mm < 0 ){
+    haha.Z();
+    double delta = length/4;
+    int sign = z > 0 ? 1 : -1;
+    z = TMath::Floor( (abs(z)-gap/2)/delta )*delta + 0.5 * delta + gap/2;
+    if( z >= 107.375 ) z = 88.625;
+    z = sign * z;
+  }
+
+  haha.SetY(y);
+  haha.SetZ(z);
+  haha.RotateZ(phi);
+
+  return haha;
+
+}
+
+
+inline void SX3::CalSX3Pos(unsigned short ID, unsigned short chUp, unsigned short chDown, unsigned short chBack, float eUp, float eDown){
+
+  hitPos.Clear();
+
+
+  if( (chUp - chDown) != 1 || (chDown % 2) != 0) return ;
+
+  int reducedID = ID % numDet;
+
+  TVector3 sa, sb;
+
+  if( ID < numDet ){ //down
+    sa = SDn[reducedID].second;
+    sb = SDn[reducedID].first;
+  }else{
+    sa = SUp[reducedID].second;
+    sb = SUp[reducedID].first;
+  }
+
+  hitPos.SetX( (sb.X() - sa.X()) * chUp/8 + sa.X());
+  hitPos.SetY( (sb.Y() - sa.Y()) * chUp/8 + sa.Y());
+
+  if( eUp == 0 || eDown == 0 ){
+    hitPos.SetZ( sa.Z() + (2*(chBk - 7)-1) * length / 8 );
+  }else{
+    double frac = (eUp - eDown)/(eUp + eDown); // from +1 (downstream) to -1 (upstream)
+    double zPos = sa.Z() +  length * frac/2;
+    hitPos.SetZ( zPos );
+  }
+
+}
+
+#endif 

Armory/Armory/ClassTransfer.h

@@ -0,0 +1,468 @@
+#ifndef ClassTransfer_h
+#define ClassTransfer_h
+
+#include "TBenchmark.h"
+#include "TLorentzVector.h"
+#include "TVector3.h"
+#include "TMath.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TRandom.h"
+#include "TMacro.h"
+#include "TGraph.h"
+#include <vector>
+#include <fstream>
+
+#include "Isotope.h"
+
+class ReactionConfig{
+
+public:
+
+  ReactionConfig(){}
+  ~ReactionConfig(){}
+
+  int beamA, beamZ;       
+  int targetA, targetZ;     
+  int recoilLightA, recoilLightZ;
+  int recoilHeavyA, recoilHeavyZ;
+
+  float beamEnergy;         ///MeV/u
+  float beamEnergySigma;    ///beam-energy_sigma_in_MeV/u
+  float beamAngle;          ///beam-angle_in_mrad
+  float beamAngleSigma;     ///beam-emittance_in_mrad
+  float beamX;              ///x_offset_of_Beam_in_mm
+  float beamY;              ///y_offset_of_Beam_in_mm
+  int numEvents;            ///number_of_Event_being_generated
+  bool isTargetScattering;  ///isTargetScattering
+  float targetDensity;      ///target_density_in_g/cm3
+  float targetThickness;    ///targetThickness_in_cm
+  std::string beamStoppingPowerFile;         ///stopping_power_for_beam
+  std::string recoilLightStoppingPowerFile;  ///stopping_power_for_light_recoil
+  std::string recoilHeavyStoppingPowerFile;  ///stopping_power_for_heavy_recoil
+  bool isDecay;        ///isDacay
+  int heavyDecayA;     ///decayNucleus_A
+  int heavyDecayZ;     ///decayNucleus_Z
+  bool isRedo;         ///isReDo
+  std::vector<float> beamEx;        ///excitation_energy_of_A[MeV]
+
+
+  void SetReaction(int beamA, int beamZ,
+                   int targetA, int targetZ,
+                   int recoilA, int recoilZ, float beamEnergy_AMeV){
+    this->beamA = beamA;
+    this->beamZ = beamZ;
+    this->targetA = targetA;
+    this->targetZ = targetZ;
+    this->recoilLightA = recoilA;
+    this->recoilLightZ = recoilZ;
+
+    recoilHeavyA = this->beamA + this->targetA - recoilLightA;
+    recoilHeavyZ = this->beamZ + this->targetZ - recoilLightZ;
+
+  } 
+
+  void LoadReactionConfig(TMacro * macro){
+
+    if( macro == NULL ) return ;
+
+    int numLine = macro->GetListOfLines()->GetSize();
+
+    for( int i = 0; i < numLine; i ++){
+
+      std::vector<std::string> str = SplitStr(macro->GetListOfLines()->At(i)->GetName(), " ");
+
+      ///printf("%d | %s\n", i,  str[0].c_str());
+
+      if( str[0].find_first_of("#") == 0 ) break;
+
+      if( i ==  0 ) beamA           = atoi(str[0].c_str());
+      if( i ==  1 ) beamZ           = atoi(str[0].c_str());
+      if( i ==  2 ) targetA         = atoi(str[0].c_str());
+      if( i ==  3 ) targetZ         = atoi(str[0].c_str());
+      if( i ==  4 ) recoilLightA    = atoi(str[0].c_str());
+      if( i ==  5 ) recoilLightZ    = atoi(str[0].c_str());
+      if( i ==  6 ) beamEnergy      = atof(str[0].c_str());
+      if( i ==  7 ) beamEnergySigma = atof(str[0].c_str());
+      if( i ==  8 ) beamAngle       = atof(str[0].c_str());
+      if( i ==  9 ) beamAngleSigma  = atof(str[0].c_str());
+      if( i == 10 ) beamX           = atof(str[0].c_str());
+      if( i == 11 ) beamY           = atof(str[0].c_str());
+      if( i == 12 ) numEvents       = atoi(str[0].c_str());
+      if( i == 13 ) {
+          if( str[0].compare("false") == 0 ) isTargetScattering = false;
+          if( str[0].compare("true")  == 0 ) isTargetScattering = true;
+      }
+      if( i == 14 ) targetDensity   = atof(str[0].c_str());
+      if( i == 15 ) targetThickness = atof(str[0].c_str());
+      if( i == 16 ) beamStoppingPowerFile = str[0];
+      if( i == 17 ) recoilLightStoppingPowerFile = str[0];
+      if( i == 18 ) recoilHeavyStoppingPowerFile = str[0];
+      if( i == 19 ) {
+          if( str[0].compare("false") == 0 ) isDecay = false;
+          if( str[0].compare("true")  == 0 ) isDecay = true;
+      }
+      if( i == 20 ) heavyDecayA = atoi(str[0].c_str());
+      if( i == 21 ) heavyDecayZ = atoi(str[0].c_str());
+      if( i == 22 ) {
+          if( str[0].compare("false") == 0 ) isRedo = false;
+          if( str[0].compare("true" ) == 0 ) isRedo = true;
+      }
+      
+      if( i >= 23) {
+          beamEx.push_back( atof(str[0].c_str()) );
+      }
+    }
+
+      recoilHeavyA = beamA + targetA - recoilLightA;
+      recoilHeavyZ = beamZ + targetZ - recoilLightZ;
+
+  }
+
+  void PrintReactionConfig(){
+
+  printf("=====================================================\n");
+  printf("   beam : A = %3d, Z = %2d \n", beamA, beamZ);
+  printf(" target : A = %3d, Z = %2d \n", targetA, targetZ);
+  printf("  light : A = %3d, Z = %2d \n", recoilLightA, recoilLightZ);
+
+  printf(" beam Energy : %.2f +- %.2f MeV/u, dE/E = %5.2f %%\n", beamEnergy, beamEnergySigma, beamEnergySigma/beamEnergy);
+  printf("       Angle : %.2f +- %.2f mrad\n", beamAngle, beamAngleSigma);
+  printf("      offset : (x,y) = (%.2f, %.2f) mmm \n", beamX, beamY);
+
+  printf("##### number of Simulation Events : %d \n", numEvents);
+  
+  printf("    is target scattering : %s \n", isTargetScattering ? "Yes" : "No");
+
+  if(isTargetScattering){
+    printf(" target density : %.f g/cm3\n", targetDensity);
+    printf("      thickness : %.f cm\n", targetThickness);
+    printf("         beam stopping file : %s \n", beamStoppingPowerFile.c_str());
+    printf(" recoil light stopping file : %s \n", recoilLightStoppingPowerFile.c_str());
+    printf(" recoil heavy stopping file : %s \n", recoilHeavyStoppingPowerFile.c_str());
+  }
+  
+  printf("       is simulate decay : %s \n", isDecay ? "Yes" : "No");
+  if( isDecay ){
+    printf(" heavy decay : A = %d, Z = %d \n", heavyDecayA, heavyDecayZ);
+  }
+  printf(" is Redo until hit array : %s \n", isRedo ? "Yes" : "No");
+
+  printf(" beam Ex : %.2f MeV \n", beamEx[0]);
+  for( int i = 1; i < (int) beamEx.size(); i++){
+    printf("         %.2f MeV \n", beamEx[i]);
+  }
+  
+  printf("=====================================================\n");
+
+  }
+
+};
+
+//=======================================================
+//#######################################################
+// Class for Transfer Reaction
+// reaction notation A(a,b)B
+// A = incident particle
+// a = target
+// b = light scattered particle
+// B = heavy scattered particle
+//======================================================= 
+class TransferReaction {
+public:
+  TransferReaction();
+  ~TransferReaction();
+
+  void SetA(int A, int Z, double Ex);
+  void Seta(int A, int Z);
+  void Setb(int A, int Z);
+  void SetB(int A, int Z);
+  void SetIncidentEnergyAngle(double KEA, double theta, double phi);
+  void SetExA(double Ex);
+  void SetExB(double Ex);
+  void SetReactionFromFile(string settingFile);
+  
+  TString GetReactionName();
+  TString GetReactionName_Latex();
+  
+  ReactionConfig GetRectionConfig() { return reaction;}
+  
+  double GetMass_A(){return mA + ExA;}
+  double GetMass_a(){return ma;}
+  double GetMass_b(){return mb;}
+  double GetMass_B(){return mB + ExB;}
+  
+  double GetCMTotalKE() {return Etot - mA - ma;}
+  double GetQValue()    {return mA + ExA + ma - mb - mB - ExB;}
+  double GetMaxExB()    {return Etot - mb - mB;}
+  
+  TLorentzVector GetPA(){return PA;}
+  TLorentzVector GetPa(){return Pa;}
+  TLorentzVector GetPb(){return Pb;}
+  TLorentzVector GetPB(){return PB;}
+  
+  void CalReactionConstant();
+  
+  TLorentzVector * Event(double thetaCM, double phiCM);
+  
+  double GetEx(){return Ex;}
+  double GetThetaCM(){return thetaCM;}
+  
+  double GetMomentumbCM()   {return p;}
+  double GetReactionBeta()  {return beta;}
+  double GetReactionGamma() {return gamma;}
+  double GetCMTotalEnergy() {return Etot;}
+  
+private:
+
+  ReactionConfig reaction;
+
+  string nameA, namea, nameb, nameB;
+  double thetaIN, phiIN;
+  double mA, ma, mb, mB;
+
+  double TA, T; // TA = KE of A pre u, T = total energy
+  double ExA, ExB;
+  double Ex, thetaCM; //calculated Ex using inverse mapping from e and z to thetaCM
+  
+  bool isReady;
+  bool isBSet;
+  
+  double k; // CM Boost momentum
+  double beta, gamma; //CM boost beta
+  double Etot;
+  double p; // CM frame momentum of b, B
+  
+  TLorentzVector PA, Pa, Pb, PB;
+
+  TString format(TString name);
+   
+};
+
+TransferReaction::TransferReaction(){
+  
+   thetaIN = 0.;
+   phiIN = 0.;
+   SetA(24, 12, 0);
+   Seta(4,2);
+   Setb(1,1);
+   SetB(27,13);
+   TA = 2.5;
+   T = TA * reaction.beamA;
+   
+   ExA = 0;
+   ExB = 0;
+   
+   Ex = TMath::QuietNaN();
+   thetaCM = TMath::QuietNaN();
+   
+   CalReactionConstant();
+   
+   TLorentzVector temp (0,0,0,0);
+   PA = temp;
+   Pa = temp;
+   Pb = temp;
+   PB = temp;
+   
+}
+
+TransferReaction::~TransferReaction(){
+
+}
+
+void TransferReaction::SetA(int A, int Z, double Ex = 0){
+  Isotope temp (A, Z);
+  mA = temp.Mass;
+  reaction.beamA = A;
+  reaction.beamZ = Z;
+  ExA = Ex;
+  nameA = temp.Name;
+  isReady = false;
+  isBSet = true;
+  
+}
+
+void TransferReaction::Seta(int A, int Z){
+  Isotope temp (A, Z);
+  ma = temp.Mass;
+  reaction.targetA = A;
+  reaction.targetZ = Z;
+  namea = temp.Name;
+  isReady = false;
+  isBSet = false;
+}
+
+void TransferReaction::Setb(int A, int Z){
+  Isotope temp (A, Z);
+  mb = temp.Mass;
+  reaction.recoilLightA = A;
+  reaction.recoilLightZ = Z;
+  nameb = temp.Name;
+  isReady = false;
+  isBSet = false;
+}
+void TransferReaction::SetB(int A, int Z){
+  Isotope temp (A, Z);
+  mB = temp.Mass;
+  reaction.recoilHeavyA = A;
+  reaction.recoilHeavyZ = Z;
+  nameB = temp.Name;
+  isReady = false;
+  isBSet = true;
+}
+
+void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double phi){
+  this->TA = KEA;
+  this->T = TA * reaction.beamA;
+  this->thetaIN = theta;
+  this->phiIN = phi;
+  isReady = false;
+}
+
+void TransferReaction::SetExA(double Ex){
+  this->ExA = Ex;
+  isReady = false;
+}
+
+void TransferReaction::SetExB(double Ex){
+  this->ExB = Ex;
+  isReady = false;
+}
+
+void TransferReaction::SetReactionFromFile(string settingFile){
+
+   TMacro * haha = new TMacro();
+   if( haha->ReadFile(settingFile.c_str()) > 0 ) {
+      reaction.LoadReactionConfig(haha);
+
+     SetA(reaction.beamA, reaction.beamZ);
+     Seta(reaction.targetA, reaction.targetZ);
+     Setb(reaction.recoilLightA, reaction.recoilLightZ);
+     SetB(reaction.recoilHeavyA, reaction.recoilHeavyZ);
+
+     SetIncidentEnergyAngle(reaction.beamEnergy, 0, 0);
+     CalReactionConstant();
+   }else{
+     
+     printf("cannot read file %s.\n", settingFile.c_str());
+     isReady = false;
+   }
+   
+}
+
+TString TransferReaction::GetReactionName(){
+  TString rName;
+  rName.Form("%s(%s,%s)%s", nameA.c_str(), namea.c_str(), nameb.c_str(), nameB.c_str()); 
+  return rName;
+}
+
+TString TransferReaction::format(TString name){
+  if( name.IsAlpha() ) return name;
+  int len = name.Length();
+  TString temp = name;
+  TString temp2 = name;
+  if( temp.Remove(0, len-2).IsAlpha()){
+     temp2.Remove(len-2);
+  }else{
+     temp = name;
+     temp.Remove(0, len-1);
+     temp2.Remove(len-1);
+  }
+  return "^{"+temp2+"}"+temp;
+}
+
+TString TransferReaction::GetReactionName_Latex(){
+  TString rName;
+  rName.Form("%s(%s,%s)%s", format(nameA).Data(), format(namea).Data(), format(nameb).Data(), format(nameB).Data()); 
+  return rName;
+}
+
+void TransferReaction::CalReactionConstant(){
+   if( !isBSet){
+      reaction.recoilHeavyA = reaction.beamA + reaction.targetA - reaction.recoilLightA;
+      reaction.recoilHeavyZ = reaction.beamZ + reaction.targetZ - reaction.recoilLightZ;
+      Isotope temp (reaction.recoilHeavyA, reaction.recoilHeavyZ);
+      mB = temp.Mass;
+      isBSet = true;
+   }
+   
+   k = TMath::Sqrt(TMath::Power(mA + ExA + T, 2) - (mA + ExA) * (mA + ExA)); 
+   beta = k / (mA + ExA + ma + T);
+   gamma = 1 / TMath::Sqrt(1- beta * beta);   
+   Etot = TMath::Sqrt(TMath::Power(mA + ExA + ma + T,2) - k * k);
+   p = TMath::Sqrt( (Etot*Etot - TMath::Power(mb + mB + ExB,2)) * (Etot*Etot - TMath::Power(mb - mB - ExB,2)) ) / 2 / Etot;
+   
+   PA.SetXYZM(0, 0, k, mA + ExA);
+   PA.RotateY(thetaIN);
+   PA.RotateZ(phiIN);
+   
+   Pa.SetXYZM(0,0,0,ma);
+   
+   isReady = true;
+}
+
+TLorentzVector * TransferReaction::Event(double thetaCM, double phiCM)
+{
+   if( isReady == false ){
+      CalReactionConstant();
+   }
+
+   //TLorentzVector Pa(0, 0, 0, ma);
+   
+   //---- to CM frame
+   TLorentzVector Pc = PA + Pa;
+   TVector3 b = Pc.BoostVector();
+   
+   TVector3 vb(0,0,0);
+   
+   if( b.Mag() > 0 ){
+     TVector3 v0 (0,0,0);
+     TVector3 nb = v0 - b;
+     
+     TLorentzVector PAc = PA; 
+     PAc.Boost(nb);
+     TVector3 vA = PAc.Vect();
+     
+     TLorentzVector Pac = Pa;
+     Pac.Boost(nb);
+     TVector3 va = Pac.Vect();
+     
+     //--- construct vb
+     vb = va;
+     vb.SetMag(p);
+
+     TVector3 ub = vb.Orthogonal();
+     vb.Rotate(thetaCM, ub);
+     vb.Rotate(phiCM + TMath::PiOver2(), va); // somehow, the calculation turn the vector 90 degree.
+     //vb.Rotate(phiCM , va); // somehow, the calculation turn the vector 90 degree.
+   }
+   
+   //--- from Pb
+   TLorentzVector Pbc;
+   Pbc.SetVectM(vb, mb);
+   
+   //--- from PB
+   TLorentzVector PBc;
+   //PBc.SetVectM(vB, mB + ExB);
+   PBc.SetVectM(-vb, mB + ExB);
+   
+   //---- to Lab Frame
+   TLorentzVector Pb = Pbc;
+   Pb.Boost(b);
+   TLorentzVector PB = PBc;
+   PB.Boost(b);
+   
+   TLorentzVector * output = new TLorentzVector[4];
+   output[0] = PA;
+   output[1] = Pa;
+   output[2] = Pb;
+   output[3] = PB;
+   
+   this->Pb = Pb;
+   this->PB = PB;
+   
+   
+   return output;   
+}
+
+#endif 

Armory/Armory/HistPlotter.h

@@ -0,0 +1,414 @@
+#ifndef HISTPLOTTER_H
+#define HISTPLOTTER_H
+#include <TCanvas.h>
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TStyle.h>
+#include <iostream>
+#include <TFile.h>
+#include <TMemFile.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TCutG.h>
+#include <signal.h>
+#include <cstdlib>
+#include <utility>
+#include <fstream>
+#include <sstream>
+#include <unordered_map>
+#include <set>
+#include <TGraphErrors.h>
+
+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::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
+    TFile *ofile=nullptr; //!< TFile pointer for the output file
+    TMemFile *omfile=nullptr; //!< TFile pointer for the output memfile
+
+    //Caches to permit FillN() calls
+    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);
+    inline void FlushToDisk(); //!< Writes all objects to file before closing, nesting objects in folders as is found necessary
+    inline void PrintObjects(); //!< Dump objects to std::cout for inspection
+    inline void ReadCuts(std::string); 
+    inline TCutG* FindCut(std::string cut) {
+        return static_cast<TCutG*>(cutsMap.at(cut));
+    }
+    inline void set_barrier_limit(long long limit) { barrier_limit = limit; }
+    inline void barrier_increment() {
+        barrier_count++;
+        if(barrier_count == barrier_limit) {
+            FillN_All_Histograms();
+            barrier_count=0;
+        }
+    }
+    /*! \fn void FindCut()
+        \brief        
+            - Searches for a cut by name 'cut' in the internal list of cuts 'cutsMap'. Ugly fails (via unresolved at()) if such a cut isn't found.            
+        \param filename - name of the plainxtext file containing the cut file locations and identifiers
+        \return Pointer to the TCutG object that matches the name. Very useful to use this as plotter.FindCut("protonbarrelpid")->IsInside(deltaE, E) for instance.
+    */
+
+    inline void SetNewTitle(std::string name, std::string title) {
+        auto result = oMap.find(name); //result is an iterator
+        if(result==oMap.end()) return; //no warnings, could be changed in future
+        else
+        	static_cast<TNamed*>(oMap.at(name))->SetTitle(title.c_str()); // set new title
+	}
+
+    //Smart functions that create a new histogram if it doesn't exist.
+    inline void FillGraph(const std::string &name, float valuex, float valuey, float errx=0, float erry=0);
+    inline void Fill1D(const std::string& name,int nbinsx, float xlow, float xhigh, float value);
+    inline void Fill2D(const std::string& name,int nbinsx, float xlow, float xhigh
+                ,int nbinsy, float ylow, float yhigh, float valuex, float valuey);
+    inline void Fill1D(const std::string& name,int nbinsx, float xlow, float xhigh, float value, const std::string& folder);
+    inline void Fill2D(const std::string& name,int nbinsx, float xlow, float xhigh
+                ,int nbinsy, float ylow, float yhigh, float valuex, float valuey, const std::string& folder);
+    //TObject* findObject(std::string key);
+};
+
+HistPlotter::HistPlotter(std::string outfile, std::string type="") {
+    /*!
+        \brief Constructor. Opens a TFile instance with the specified filename
+        \param outfile : std::string that holds the desired output ROOT filename
+        \return None
+    */
+    if(type=="" || type == "TFILE") {
+        ofile = new TFile(outfile.c_str(),"recreate");
+        filetype = TFILE;
+    } else if(type =="TMEMFILE") {
+        omfile = new TMemFile(outfile.c_str(),"recreate");
+        filetype=TMEMFILE;
+    } else {
+        std::cout << "Unknown type "<< type << " specified for HistPlotter (use \"TFILE\" or \"TMEMFILE\"), using default \"TFILE\"  " << std::endl;
+        ofile = new TFile(outfile.c_str(),"recreate");
+        filetype = TFILE;
+    }
+    barrier_count=0;
+    barrier_limit=1000;
+}
+
+void HistPlotter::FillN_All_Histograms() {
+    for(auto it=oMap.begin(); it!=oMap.end(); it++ ) {
+        //it->first is std::string 'name', it->second is the TObject
+        if(it->second->InheritsFrom("TH1F")) {
+            //FillN(size, array-of-doubles, array-of-weights); //we set array-of-weights to (1,1,1,.. (size)
+            static_cast<TH1F*>(it->second)->FillN(onedimcache[it->first].size(), //size
+                                                  onedimcache[it->first].data(), //array
+                                                  std::vector<double>(onedimcache[it->first].size(),1.0).data()); //weight of ones
+            onedimcache[it->first].clear();
+        }  else if(it->second->InheritsFrom("TH2F")) {
+            //FillN(size, array-of-doubles, array-of-weights); //we set array-of-weights to (1,1,1,.. (size))
+            static_cast<TH2F*>(it->second)->FillN(twodimcache[it->first].first.size(), //size
+                                                  twodimcache[it->first].first.data(), //x array
+                                                  twodimcache[it->first].second.data(), //y array
+                                                  std::vector<double>(twodimcache[it->first].first.size(),1.0).data()); //weight of ones
+           twodimcache[it->first].first.clear();
+           twodimcache[it->first].second.clear();
+        }
+    }
+    std::cout << "." << std::endl;
+}
+void HistPlotter::FlushToDisk() {
+    /*! \fn void FlushToDisk()
+        \brief         Function that can be used at any point to exit smoothly by saving all ROOT objects in memory
+        to the output file before closing it. Obeys the binding of histograms to separate folders, if so specified.
+        \return No return -- void
+    */
+    if(filetype==TMEMFILE && omfile) {
+        std::cout << "Not flushing a TMemfile .. exiting .." << std::endl;
+        delete omfile;
+        return;
+    }
+    if(ofile->IsZombie() || !ofile) {
+        std::cerr << "Output file is zombie, finishing up without writing to disk!" << std::endl;
+        return;
+    }
+    FillN_All_Histograms();
+    for(auto it=oMap.begin(); it!=oMap.end(); it++ ) {
+        //omap maps: name(first) to object address(second).
+        // foldersForObjects maps: object address(first) to foldername(second)
+        auto result = foldersForObjects.find(it->second); //returns <TObject* histogram,std::string foldername> pair if found
+        if(result!=foldersForObjects.end()) { //we try to create folder if needed and cd to it
+            ofile->mkdir(result->second.c_str(),"",kTRUE); // args: name, title, returnExistingDirectory
+            ofile->cd(result->second.c_str());
+        } else {
+            ofile->cd(); //toplevel for all default histograms. Default setting
+        }
+        it->second->Write();
+    }
+
+    //Create a directory for all cuts, and save all cuts in them
+    ofile->mkdir("gCUTS","",kTRUE);
+    ofile->cd("gCUTS");
+    for(auto it=cutsMap.begin(); it!=cutsMap.end(); it++) {
+        (static_cast<TNamed*>(it->second))->SetName(it->first.c_str());
+        it->second->Write();
+    }
+    ofile->Close();
+    std::cout << "Wrote " << oMap.size() << " histograms to TFile " << std::string(ofile->GetName()) << std::endl;
+}
+
+void HistPlotter::FillGraph(const std::string& name, float valuex, float valuey, float errx, float erry) {
+    /*! \fn void FillGraph()
+        \brief
+        - Creates a TGraphError in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
+        - Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TGraph*
+
+        \param name name of the TGraph
+        \param valuex The xvalue
+        \param valuey The yvalue
+        \param errx The x error
+        \param erry The y error
+        \return No return void
+    */
+    auto result = oMap.find(name); 
+    if(result==oMap.end()) {
+        TGraphErrors *tempG = new TGraphErrors();
+        tempG->SetName(name.c_str());
+        oMap.insert(std::make_pair(name,static_cast<TObject*>(tempG)));
+    }
+    if(!oMap.at(name)->InheritsFrom("TGraphErrors")) {
+        std::cerr << "Object " << name << " refers to something other than a TGraph*, not filling it hence!" << std::endl;
+        std::cerr << "Abort.." << std::endl;
+        FlushToDisk();
+        exit(-1);
+    }
+    // static_cast<TGraphErrors*>(oMap.at(name))->AddPointError(valuex,valuey,errx,erry);
+}
+
+void HistPlotter::Fill1D(const std::string& name, int nbinsx, float xlow, float xhigh, float value) {
+    /*! \fn void Fill1D()
+        \brief  
+        - Creates a TH1F in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
+        - Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TH1*
+
+        \param name name of the TH1F histogram
+        \param nbinsx Number of bins in the histogram
+        \param xlow Lower limit on x-axis
+        \param xhigh Upper limit on x-axis
+        \param value The bin corresponding to value in (nbinsx, xlow, xhigh) is incremented by 1
+        \return No return void
+    */
+    auto result = oMap.find(name); //result is an iterator
+    if(result==oMap.end()) {
+        TH1F* temp1D = new TH1F(name.c_str(), name.c_str(), nbinsx, xlow, xhigh);
+        oMap.insert(std::make_pair(name,static_cast<TObject*>(temp1D)));
+        onedimcache.insert(std::make_pair(name, std::vector<double>()));
+        onedimcache[name].reserve(16384);
+    } else if(foldersForObjects.find(oMap.at(name))!=foldersForObjects.end()) { //shouldn't have a folder associated with it
+            std::cerr << "Object " << name << " already registered at " << foldersForObjects[oMap[name]] <<  ", choose a different name for the histogram to be stored in toplevel .." << std::endl;
+    }
+
+    //Check if the string 'name' maps to a 1D hist. If there's any other object by this name raise issue
+    if(!oMap.at(name)->InheritsFrom("TH1F")) {
+        std::cerr << "Object " << name << " refers to something other than a TH1*, not filling it hence!" << std::endl;
+        std::cerr << "Abort.." << std::endl;
+        FlushToDisk();
+        exit(-1);
+    }
+    onedimcache[name].emplace_back(value);
+    //static_cast<TH1F*>(oMap.at(name))->Fill(value);
+}
+
+void HistPlotter::Fill1D(const std::string& name, int nbinsx, float xlow, float xhigh, float value, const std::string& foldername) {
+    /*! \fn void Fill1D()
+        \brief        
+      -  Creates a TH1F in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
+      -  Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TH1*
+      -  Remembers the foldername this particular histogram maps to, if provided. If not, defaults to toplevel.
+
+        \param name name of the TH1F histogram
+        \param nbinsx Number of bins in the histogram
+        \param xlow Lower limit on x-axis
+        \param xhigh Upper limit on x-axis
+        \param value The bin corresponding to value in (nbinsx, xlow, xhigh) is incremented by 1
+        \param foldername Name of the folder to put this histogram into. Defaults to toplevel if left empty
+        \return No return -- void
+    */
+
+    auto result = oMap.find(name); //result is an iterator
+    if(result==oMap.end()) {
+        TH1F* temp1D = new TH1F(name.c_str(), name.c_str(), nbinsx, xlow, xhigh);
+        oMap.insert(std::make_pair(name,static_cast<TObject*>(temp1D)));
+        onedimcache.insert(std::make_pair(name, std::vector<double>()));
+        onedimcache[name].reserve(16384);
+        if(foldername!="") {
+            if(folderList.find(foldername)==folderList.end()) {
+                folderList.insert(foldername);
+            }
+            foldersForObjects.insert(std::make_pair(static_cast<TObject*>(temp1D),foldername));
+        }
+    } else {
+        //object is present in map, but we enforce unique names
+        //it must already have a folder attached to it
+        if(foldersForObjects.find(oMap.at(name))==foldersForObjects.end()) {
+            std::cerr << "Object " << name << " already registered at toplevel, choose a different name for the histogram to be stored in " << foldername << " folder .." << std::endl;
+        } else if(foldersForObjects[oMap[name]]!=foldername) {
+            std::cerr << "Object " << name << " already registered at " << foldersForObjects[oMap[name]] <<  ", choose a different name for the histogram to be stored in " << foldername << " folder .." << std::endl;
+        }
+    }
+    //Check if the string 'name' maps to a 1D hist. If there's any other object by this name raise issue
+    if(!oMap.at(name)->InheritsFrom("TH1F")) {
+        std::cerr << "Object " << name << " refers to something other than a TH1*, not filling it hence!" << std::endl;
+        std::cerr << "Abort.." << std::endl;
+        FlushToDisk();
+        exit(-1);
+    }
+    onedimcache[name].emplace_back(value);
+    //static_cast<TH1F*>(oMap.at(name))->Fill(value);
+}
+
+void HistPlotter::Fill2D(const std::string& name, int nbinsx, float xlow, float xhigh, int nbinsy, float ylow, float yhigh, float valuex, float valuey) {
+    /*! \fn void Fill2D()
+        \brief        
+      -  Creates a TH2F in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
+      -  Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TH2*
+        \param name name of the TH1F histogram
+        \param nbinsx Number of xbins in the histogram
+        \param xlow Lower limit on x-axis
+        \param xhigh Upper limit on x-axis
+        \param nbinsy Number of ybins in the histogram
+        \param ylow Lower limit on y-axis
+        \param yhigh Upper limit on y-axis
+        \param valuex
+        \param valuey The bin corresponding to (valuex, valuey) in (nbinsx, xlow, xhigh, ybinsx, ylow, yhigh) is incremented by 1
+        \return No return -- void
+    */
+
+    auto result = oMap.find(name); //result is an iterator
+    if(result==oMap.end()) {
+        TH2F* temp2D = new TH2F(name.c_str(), name.c_str(), nbinsx, xlow, xhigh, nbinsy, ylow, yhigh);
+        oMap.insert(std::make_pair(name,static_cast<TObject*>(temp2D)));
+        twodimcache.insert(std::make_pair(name, std::make_pair(std::vector<double>(),std::vector<double>())));
+        twodimcache[name].first.reserve(16384);
+        twodimcache[name].second.reserve(16384);
+    } else if(foldersForObjects.find(oMap.at(name))!=foldersForObjects.end()) { //shouldn't have a folder associated with it
+        std::cerr << "Object " << name << " already registered at " << foldersForObjects[oMap[name]] <<  ", choose a different name for the histogram to be stored in toplevel .." << std::endl;
+    }
+
+    //Check if the string 'name' maps to a 1D hist. If there's any other object by this name raise issue
+    if(!oMap.at(name)->InheritsFrom("TH2F")) {
+        std::cerr << "Object " << name << " refers to something other than a TH2*, not filling it hence!" << std::endl;
+        std::cerr << "Abort.." << std::endl;
+        FlushToDisk();
+        exit(-1);
+    }
+    twodimcache[name].first.emplace_back(valuex);
+    twodimcache[name].second.emplace_back(valuey);
+    //static_cast<TH2F*>(oMap.at(name))->Fill(valuex,valuey);
+}
+
+void HistPlotter::Fill2D(const std::string& name, int nbinsx, float xlow, float xhigh, int nbinsy, float ylow, float yhigh, float valuex, float valuey, const std::string& foldername) {
+    /*! \fn void Fill2D()
+        \brief        
+      -  Creates a TH2F in memory with name 'name' if it doesn't exist, and fills it with valuex, valuey
+      -  Writes present state to disk and fails with return value -1 if the name clashes with another object that's not of type TH2*
+      -  Remembers the foldername this particular histogram maps to, if provided. If not defaults to toplevel
+
+        \param name name of the TH1F histogram
+        \param nbinsx Number of xbins in the histogram
+        \param xlow Lower limit on x-axis
+        \param xhigh Upper limit on x-axis
+        \param nbinsy Number of ybins in the histogram
+        \param ylow Lower limit on y-axis
+        \param yhigh Upper limit on y-axis
+        \param valuex
+        \param valuey The bin corresponding to (valuex, valuey) in (nbinsx, xlow, xhigh, ybinsx, ylow, yhigh) is incremented by 1
+        \param foldername Name of the folder to put this histogram into. Defaults to toplevel if left empty
+        \return No return -- void
+    */
+
+    auto result = oMap.find(name); //result is an iterator
+    if(result==oMap.end()) {
+        TH2F* temp2D = new TH2F(name.c_str(), name.c_str(), nbinsx, xlow, xhigh, nbinsy, ylow, yhigh);
+        oMap.insert(std::make_pair(name,static_cast<TObject*>(temp2D)));
+        twodimcache.insert(std::make_pair(name, std::make_pair(std::vector<double>(),std::vector<double>())));
+        twodimcache[name].first.reserve(16384);
+        twodimcache[name].second.reserve(16384);
+        if(foldername!="") {
+            if(folderList.find(foldername)==folderList.end()) {
+                folderList.insert(foldername);
+            }
+            foldersForObjects.insert(std::make_pair(static_cast<TObject*>(temp2D),foldername));
+        }
+    } else {
+        //object is present in map, but we enforce unique names
+        //it must already have a folder attached to it
+        if(foldersForObjects.find(oMap.at(name))==foldersForObjects.end()) {
+            std::cerr << "Object " << name << " already registered at toplevel, choose a different name for the histogram to be stored in " << foldername << " folder .." << std::endl;
+        } else if(foldersForObjects[oMap.at(name)]!=foldername) {
+            std::cerr << "Object " << name << " already registered at " << foldersForObjects[oMap[name]] <<  ", choose a different name for the histogram to be stored in " << foldername << " folder .." << std::endl;
+        }
+    }
+
+    //Check if the string 'name' maps to a 1D hist. If there's any other object by this name raise issue
+    if(!oMap.at(name)->InheritsFrom("TH2F")) {
+        std::cerr << "Object " << name << " refers to something other than a TH2*, not filling it hence!" << std::endl;
+        std::cerr << "Abort.." << std::endl;
+        FlushToDisk();
+        exit(-1);
+    }
+    twodimcache[name].first.emplace_back(valuex);
+    twodimcache[name].second.emplace_back(valuey);
+    //static_cast<TH2F*>(oMap.at(name))->Fill(valuex,valuey);
+}
+
+void HistPlotter::ReadCuts(std::string filename) {
+    /*! \fn void ReadCuts()
+        \brief  Reads a list of cuts from a file. The file must have the format below, two columns
+                   - Column#1 - path to a file that contains a single TCutG object named "CUTG", the default name in ROOT.
+                   - Column#2 - The identifier name you plan to use in the code, like 'protonbarrelpid' or something, that will be searched by FindCut()
+        \param filename name of the plainxtext file containing the cut file locations and identifiers
+        \return No return -- void
+    */
+
+    std::ifstream infile;
+    infile.open(filename);
+    std::string cutfilename, cutname;
+    for(std::string line; std::getline(infile, line); ) {
+        if(line.size()!=0 && line[0]=='#')
+            ; //don't do anything with '#' lines
+        else {
+            std::stringstream ss(line);
+            ss>>cutfilename>>cutname;
+
+            TFile f(cutfilename.c_str());
+            if(f.IsZombie()) {
+                std::cerr << "Cannot open cutfile " << cutfilename << " .. skipping.." << std::endl;
+                continue;
+            }
+            TCutG *cut = (TCutG*)(f.Get("CUTG"));
+            cutsMap.insert(std::make_pair(cutname,static_cast<TObject*>(cut)));
+            f.Close();
+        } //else
+    }//for loop
+    infile.close();
+}
+
+void HistPlotter::PrintObjects() {
+    /*
+        void PrintObjects()
+        Prints the contents of the unordered_maps oMap and cutsMap to facilitate debugging
+
+    */
+    std::cout << "Type | Name " << std::endl;
+    std::cout << "---- | --------------------- " << std::endl;
+    for(auto it=oMap.begin(); it!=oMap.end(); it++ ) {
+        std::cout << it->second->ClassName() << " | "<< it->first << std::endl;
+    }
+    for(auto it=cutsMap.begin(); it!=cutsMap.end(); it++ ) {
+        std::cout << it->second->ClassName() << " | "<< it->first << std::endl;
+    }
+    std::cout << "---- | --------------------- " << std::endl;
+}
+
+#endif

Armory/Armory/Isotope.h

@@ -0,0 +1,522 @@
+/***********************************************************************
+ * 
+ *  This is Isotope.h, To extract the isotope mass from massXX.txt
+ * 
+ *-------------------------------------------------------
+ *  created by Ryan (Tsz Leung) Tang, Nov-18, 2018
+ *  email: goluckyryan@gmail.com
+ * ********************************************************************/
+
+
+#ifndef ISOTOPE_H
+#define ISOTOPE_H
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <stdio.h>
+#include <stdlib.h>
+#include "constant.h" // amu
+#include <stdlib.h>  //atoi
+#include <algorithm>
+using namespace std;
+
+string massData="mass20.txt";
+
+// about the mass**.txt
+// Mass Excess = (ATOMIC MASS - A)*amu | e.g. n : (1.088664.91585E-6-1)*amu
+// mass excess uncertaintly
+// BEA = (Z*M(1H) + N*M(1n) - Me(A,Z))/A , Me is the mass with electrons
+// BEA = (Z*mp + N*mn - M(A,Z))/A , M is the mass without electrons
+
+class Isotope {
+public:
+  int A, Z;
+  double Mass, MassError, BEA;
+  string Name, Symbol;
+  string dataSource;
+  
+  Isotope(){ dataSource = massData; };
+  Isotope(int a, int z){ dataSource = massData; SetIso(a,z);  };
+  Isotope(string name){ dataSource = massData; SetIsoByName(name); };
+
+  void SetIso(int a, int z);
+  void SetIsoByName(string name);
+
+  double CalSp(int Np, int Nn); // this for the Np-proton, Nn-neutron removal 
+  double CalSp2(int a, int z); // this is for (a,z) nucleus removal
+
+  double CalBeta(double T){
+    // double Etot = Mass + T;
+    double gamma = 1 + T/Mass;
+    double beta = sqrt(1 - 1 / gamma / gamma ) ;
+    return beta;
+  }
+
+  void Print();
+  void ListShell();
+   
+private:
+  void FindMassByAZ(int a, int z); // give mass, massError, BEA, Name, Symbol;
+  void FindMassByName(string name); // give Z, mass, massError, BEA;
+
+  int TwoJ(int nShell);
+  string Orbital(int nShell);
+  int magic(int i){
+    switch (i){
+      case 0: return   2; break;
+      case 1: return   8; break;
+      case 2: return  20; break;
+      case 3: return  28; break;
+      case 4: return  40; break;
+      case 5: return  50; break;
+      case 6: return  82; break;
+      case 7: return 128; break;
+    }
+    return 0;
+  }
+
+  int magicShellID(int i){
+    switch (i){
+      case 0: return   0; break;
+      case 1: return   2; break;
+      case 2: return   5; break;
+      case 3: return   6; break;
+      case 4: return   9; break;
+      case 5: return  10; break;
+      case 6: return  15; break;
+      case 7: return  21; break;
+    }
+    return 0;
+  }
+  
+  int fileStartLine;
+  int fileEndLine;
+  int lineMass050_099;
+  int lineMass100_149;
+  int lineMass150_199;
+  int lineMass200;
+  
+  
+  void setFileLines(){
+    fileStartLine = 37;
+    fileEndLine = 3594;
+    
+    lineMass050_099 = 466;
+    lineMass100_149 = 1160;
+    lineMass150_199 = 1994;
+    lineMass200     = 2774;
+  }
+  
+  char * heliosPath;
+  bool isFindOnce;
+  
+  
+};
+
+inline void Isotope::SetIso(int a, int z){
+    this->A = a;
+    this->Z = z;
+    FindMassByAZ(a,z); 
+}
+
+inline void Isotope::SetIsoByName(string name){
+    FindMassByName(name); 
+}
+
+inline void Isotope::FindMassByAZ(int A, int Z){
+  string line;
+  int    lineNum=0;
+  int    list_A, list_Z;
+
+  ifstream myfile;
+  int    flag=0;
+
+  setFileLines();
+
+  int numLineStart = fileStartLine;
+  int numLineEnd  = fileEndLine;
+
+  if ( A >= 50 && A < 100) numLineStart = lineMass050_099;
+  if ( A >=100 && A < 150) numLineStart = lineMass100_149;
+  if ( A >=150 && A < 200) numLineStart = lineMass150_199;
+  if ( A >=200 ) numLineStart           = lineMass200;
+
+  myfile.open(dataSource.c_str());
+
+  if (myfile.is_open()) {
+    while (/*! myfile.eof() &&*/ flag == 0 && lineNum <numLineEnd){
+      lineNum ++ ;
+      //printf("%3d  ",lineNum);
+      getline (myfile,line);
+
+      if (lineNum >= numLineStart ){
+        list_Z     = atoi((line.substr(10,5)).c_str());
+      	list_A     = atoi((line.substr(15,5)).c_str());
+
+      	if ( A == list_A && Z == list_Z) {
+            this->BEA       = atof((line.substr(54,11)).c_str());
+      		this->Mass      = list_Z*mp + (list_A-list_Z)*mn - this->BEA/1000*list_A;
+            this->MassError = atof((line.substr(65,7)).c_str());
+            string str = line.substr(20,2);
+            str.erase(remove(str.begin(), str.end(), ' '), str.end());
+            this->Symbol    = str;
+            
+            ostringstream ss;
+            ss << A << this->Symbol;
+            this->Name      = ss.str();
+     		   flag = 1;
+      	}else if ( list_A > A) {
+          this->BEA  = -404;
+          this->Mass = -404;
+          this->MassError = -404;
+          this->Symbol = "non";
+          this->Name   = "non";
+          break;
+        }
+
+      }
+    }
+    
+    if( this->Name == "1H" ) this->Name = "p";
+    if( this->Name == "2H" ) this->Name = "d";
+    if( this->Name == "3H" ) this->Name = "t";
+    if( this->Name == "4He" ) this->Name = "a";
+    
+    myfile.close();
+  }else {
+    printf("Unable to open %s\n", dataSource.c_str());
+  }
+}
+
+inline void Isotope::FindMassByName(string name){
+
+    // done seperate the Mass number and the name 
+  if( name == "n" ) {
+    this->Name = "1n";
+    this->BEA       = 0;
+    this->Mass      = mn;
+    this->MassError = 0;
+    this->Name      = "n";
+    this->A         = 1;
+    this->Z         = 0;
+    return;
+  }
+    if( name == "p" ) name = "1H";
+    if( name == "d" ) name = "2H";
+    if( name == "t" ) name = "3H";
+    if( name == "a" ) name = "4He";
+    
+    string temp = name;
+    int lastDigit = 0;
+
+    for(int i=0; temp[i]; i++){
+      if(temp[i] == '0') lastDigit =  i; 
+      if(temp[i] == '1') lastDigit =  i; 
+      if(temp[i] == '2') lastDigit =  i; 
+      if(temp[i] == '3') lastDigit =  i; 
+      if(temp[i] == '4') lastDigit =  i; 
+      if(temp[i] == '5') lastDigit =  i; 
+      if(temp[i] == '6') lastDigit =  i; 
+      if(temp[i] == '7') lastDigit =  i; 
+      if(temp[i] == '8') lastDigit =  i; 
+      if(temp[i] == '9') lastDigit =  i; 
+    }
+
+    this->Symbol = temp.erase(0,  lastDigit +1);
+    //check is Symbol is 2 charaters, if not, add " " at the end
+    if( this->Symbol.length() == 1 ){
+      this->Symbol = this->Symbol + " ";
+    }
+
+
+    temp = name;
+    int len = temp.length();
+    temp = temp.erase(lastDigit+1, len);
+    
+    this->A = atoi(temp.c_str());    
+    //printf(" Symbol = |%s| , Mass = %d\n", this->Symbol.c_str(), this->A);
+
+    // find the nucleus in the data
+    string line;
+    int    lineNum=0;
+    int    list_A;
+    string list_symbol;
+
+    ifstream myfile;
+    int    flag=0;
+
+    setFileLines();
+
+    int numLineStart = fileStartLine;
+    int numLineEnd  = fileEndLine;
+
+    if ( A >= 50 && A < 100) numLineStart = lineMass050_099;
+    if ( A >=100 && A < 150) numLineStart = lineMass100_149;
+    if ( A >=150 && A < 200) numLineStart = lineMass150_199;
+    if ( A >=200 ) numLineStart           = lineMass200;
+
+    myfile.open(dataSource.c_str());
+
+    if (myfile.is_open()) {
+      while (/*! myfile.eof() &&*/ flag == 0 && lineNum <numLineEnd){
+        lineNum ++ ;
+        //printf("%3d  ",lineNum);
+        getline (myfile,line);
+
+        if (lineNum >= numLineStart ){
+          list_symbol  = line.substr(20,2);
+          list_A       = atoi((line.substr(15,5)).c_str());
+
+          //printf(" A = %d, Sym = |%s| \n", list_A, list_symbol.c_str());
+
+          if ( this->A == list_A &&  this->Symbol == list_symbol) {
+            this->Z         = atoi((line.substr(10,5)).c_str());
+            this->BEA       = atof((line.substr(54,11)).c_str());
+      		this->Mass      = this->Z*mp + (list_A-this->Z)*mn - this->BEA/1000*list_A;
+            this->MassError = atof((line.substr(65,7)).c_str());
+            
+            string str = line.substr(20,2);
+            str.erase(remove(str.begin(), str.end(), ' '), str.end());
+            this->Symbol    = str;
+            
+            ostringstream ss;
+            ss << this->A << this->Symbol;
+            this->Name      = ss.str();
+            flag = 1;
+          }else if ( list_A > this->A) {
+            this->BEA  = -404;
+            this->Mass = -404;
+            this->MassError = -404;
+            this->Symbol = "non";
+            this->Name   = "non";
+            break;
+          }
+
+        }
+      }
+      myfile.close();
+    }else {
+      printf("Unable to open %s\n", dataSource.c_str());
+    }
+}
+
+inline double Isotope::CalSp(int Np, int Nn){
+  Isotope nucleusD(A - Np - Nn, Z - Np);  
+
+  if( nucleusD.Mass != -404){
+    return nucleusD.Mass + Nn*mn + Np*mp - this->Mass;
+  }else{
+    return -404;
+  }
+}
+
+inline double Isotope::CalSp2(int a, int z){
+  Isotope nucleusD(A - a , Z - z);
+  Isotope nucleusS(a,z);  
+
+  if( nucleusD.Mass != -404 && nucleusS.Mass != -404){
+    return nucleusD.Mass + nucleusS.Mass - this->Mass;
+  }else{
+    return -404;
+  }
+}
+
+inline int Isotope::TwoJ(int nShell){
+
+  switch(nShell){
+    case  0: return  1; break; // 0s1/2
+    case  1: return  3; break; // 0p3/2
+    case  2: return  1; break; // 0p1/2 -- 8
+    case  3: return  5; break; // 0d5/2
+    case  4: return  1; break; // 1s1/2
+    case  5: return  3; break; // 0d3/2 -- 20
+    case  6: return  7; break; // 0f7/2 -- 28
+    case  7: return  3; break; // 1p3/2
+    case  8: return  1; break; // 1p1/2
+    case  9: return  5; break; // 0f5/2 -- 40
+    case 10: return  9; break; // 0g9/2 -- 50
+    case 11: return  7; break; // 0g7/2
+    case 12: return  5; break; // 1d5/2
+    case 13: return 11; break; // 0h11/2
+    case 14: return  3; break; // 1d3/2
+    case 15: return  1; break; // 2s1/2 -- 82
+    case 16: return  9; break; // 0h9/2
+    case 17: return  7; break; // 1f7/2
+    case 18: return 13; break; // 0i13/2
+    case 19: return  3; break; // 2p3/2
+    case 20: return  5; break; // 1f5/2
+    case 21: return  1; break; // 1p1/2 -- 126
+    case 22: return  9; break; // 1g9/2
+    case 23: return 11; break; // 0i11/2
+    case 24: return 15; break; // 0j15/2
+    case 25: return  5; break; // 2d5/2
+    case 26: return  1; break; // 3s1/2
+    case 27: return  3; break; // 2d3/2
+    case 28: return  7; break; // 1g7/2
+  }
+
+  return 0;
+}
+
+inline string Isotope::Orbital(int nShell){
+
+  switch(nShell){
+    case  0: return  "0s1 "; break;  //
+    case  1: return  "0p3 "; break;  //
+    case  2: return  "0p1 "; break;  //-- 8
+    case  3: return  "0d5 "; break;  //
+    case  4: return  "1s1 "; break;  //
+    case  5: return  "0d3 "; break;  //-- 20
+    case  6: return  "0f7 "; break;  //-- 28
+    case  7: return  "1p3 "; break;  //
+    case  8: return  "1p1 "; break;  //
+    case  9: return  "0f5 "; break;  //-- 40
+    case 10: return  "0g9 "; break;  //-- 50
+    case 11: return  "0g7 "; break;  //
+    case 12: return  "1d5 "; break;  //
+    case 13: return  "0h11"; break;  //
+    case 14: return  "1d3 "; break;  //
+    case 15: return  "2s1 "; break;  //-- 82
+    case 16: return  "0h9 "; break;  //
+    case 17: return  "1f7 "; break;  //
+    case 18: return  "0i13"; break;  //
+    case 19: return  "2p3 "; break;  //
+    case 20: return  "1f5 "; break;  //
+    case 21: return  "1p1 "; break;  //-- 126
+    case 22: return  "1g9 "; break;  //
+    case 23: return  "0i11"; break;  //
+    case 24: return  "0j15"; break;  //
+    case 25: return  "2d5 "; break;  //
+    case 26: return  "3s1 "; break;  //
+    case 27: return  "2d3 "; break;  //
+    case 28: return  "1g7 "; break;  //
+  }
+
+  return "nan";
+}
+
+inline void Isotope::ListShell(){
+
+  if( Mass < 0 ) return;
+
+  int n = A-Z;
+  int p = Z;
+
+  int k = std::min(n,p);
+  int nMagic = 0;
+  for( int i = 0; i < 7; i++){
+    if( magic(i) < k && k <= magic(i+1) ){
+      nMagic = i;
+      break;
+    }
+  }
+
+  int coreShell = magicShellID(nMagic-1);
+  int coreZ1 = magic(nMagic-1);
+  int coreZ2 = magic(nMagic);
+
+  Isotope core1( 2*coreZ1, coreZ1);
+  Isotope core2( 2*coreZ2, coreZ2);
+
+  printf("------------------ Core:%3s, inner Core:%3s \n", (core2.Name).c_str(), (core1.Name).c_str());
+  printf("         || ");
+  int t = std::max(n,p);
+  int nShell = 0;
+  do{
+    int occ = TwoJ(nShell)+1;
+    if( nShell > coreShell) {
+      printf("%4s", Orbital(nShell).c_str());
+         if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9  || nShell == 10 || nShell == 15 || nShell == 21){
+        printf("|");
+      }else{
+        printf(",");
+      } 
+    }
+    t = t - occ;
+    nShell++;
+  }while( t > 0  && nShell < 29);
+  for( int i = 1; i <= 6; i++){
+    if (nShell < 28) {
+      printf("%4s,", Orbital(nShell).c_str());
+    }else if( nShell == 28 ) {
+      printf("%4s", Orbital(nShell).c_str());
+    }
+    nShell ++;
+  }
+  if (nShell < 29) printf("%4s", Orbital(nShell).c_str());
+  printf("\n");
+
+
+  printf(" Z = %3d || ", p);
+  nShell = 0;
+  do{
+    int occ = TwoJ(nShell)+1;
+    if( nShell > coreShell ){
+      if( p > occ ) {
+         printf("%-4d", occ);
+         if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9  || nShell == 10 || nShell == 15 || nShell == 21){
+          printf("|");
+        }else{
+          printf(",");
+        } 
+      }else{
+        printf("%-4d", p);
+      }
+    }
+    p = p - occ;
+    nShell++;
+  }while( p > 0  && nShell < 29);
+  printf("\n");
+
+  printf(" N = %3d || ", n);
+  nShell = 0;  
+  do{
+    int occ = TwoJ(nShell)+1;
+    if ( nShell > coreShell ){
+      if( n > occ ) {
+         printf("%-4d", occ);
+         if( nShell == 0 || nShell == 2 || nShell == 5 || nShell ==6 || nShell == 9  || nShell == 10 || nShell == 15 || nShell == 21){
+          printf("|");
+        }else{
+          printf(",");
+        } 
+      }else{
+        printf("%-4d", n);
+      }
+    }
+    n = n - occ;
+    nShell++;
+  }while( n > 0  && nShell < 29);
+  printf("\n");
+
+  printf("------------------ \n");
+}
+
+
+inline void Isotope::Print(){
+
+  if (Mass > 0){  
+    
+    dataSource = massData;
+      
+    printf(" using mass data : %s \n", dataSource.c_str());
+    printf(" mass of \e[47m\e[31m%s\e[m nucleus (Z,A)=(%3d,%3d) is \e[47m\e[31m%12.5f\e[m MeV, BE/A=%7.5f MeV\n", Name.c_str(), Z, A, Mass, BEA/1000.);     
+    printf(" total BE    : %12.5f MeV\n",BEA*A/1000.);
+    printf(" mass in amu : %12.5f u\n",Mass/amu);
+    printf(" mass excess : %12.5f MeV\n", Mass + Z*0.510998950 - A*amu);
+    printf("-------------- Seperation energy \n");
+    printf(" S1p: %8.4f| S1n: %8.4f| S(2H ): %8.4f| S1p1n : %8.4f\n", CalSp(1, 0), CalSp(0, 1), CalSp2(2, 1), CalSp(1, 1));
+    printf(" S2p: %8.4f| S2n: %8.4f| S(3He): %8.4f| S(3H) : %8.4f\n", CalSp(2, 0), CalSp(0, 2), CalSp2(3, 2), CalSp2(3, 1));
+    printf(" S3p: %8.4f| S3n: %8.4f| S(4He): %8.4f|\n",               CalSp(3, 0), CalSp(0, 3), CalSp2(4, 2));
+    printf(" S4p: %8.4f| S4n: %8.4f| \n",                             CalSp(4, 0), CalSp(0, 4));
+
+  }else{
+    printf("Error %6.0f, no nucleus with (Z,A) = (%3d,%3d). \n", Mass, Z, A);
+  }
+
+
+}
+
+#endif

Armory/Armory/LICENSE

@@ -0,0 +1,674 @@
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+  The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works.  By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users.  We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors.  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+  To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights.  Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received.  You must make sure that they, too, receive
+or can get the source code.  And you must show them these terms so they
+know their rights.
+
+  Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+  For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software.  For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+  Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so.  This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software.  The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable.  Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products.  If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+  Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary.  To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                       TERMS AND CONDITIONS
+
+  0. Definitions.
+
+  "This License" refers to version 3 of the GNU General Public License.
+
+  "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+  "The Program" refers to any copyrightable work licensed under this
+License.  Each licensee is addressed as "you".  "Licensees" and
+"recipients" may be individuals or organizations.
+
+  To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy.  The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+  A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+  To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy.  Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+  To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies.  Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+  An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License.  If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+  1. Source Code.
+
+  The "source code" for a work means the preferred form of the work
+for making modifications to it.  "Object code" means any non-source
+form of a work.
+
+  A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+  The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form.  A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+  The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities.  However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work.  For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+  The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+  The Corresponding Source for a work in source code form is that
+same work.
+
+  2. Basic Permissions.
+
+  All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met.  This License explicitly affirms your unlimited
+permission to run the unmodified Program.  The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work.  This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+  You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force.  You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright.  Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+  Conveying under any other circumstances is permitted solely under
+the conditions stated below.  Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+  No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+  When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+  4. Conveying Verbatim Copies.
+
+  You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+  You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+  5. Conveying Modified Source Versions.
+
+  You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+    a) The work must carry prominent notices stating that you modified
+    it, and giving a relevant date.
+
+    b) The work must carry prominent notices stating that it is
+    released under this License and any conditions added under section
+    7.  This requirement modifies the requirement in section 4 to
+    "keep intact all notices".
+
+    c) You must license the entire work, as a whole, under this
+    License to anyone who comes into possession of a copy.  This
+    License will therefore apply, along with any applicable section 7
+    additional terms, to the whole of the work, and all its parts,
+    regardless of how they are packaged.  This License gives no
+    permission to license the work in any other way, but it does not
+    invalidate such permission if you have separately received it.
+
+    d) If the work has interactive user interfaces, each must display
+    Appropriate Legal Notices; however, if the Program has interactive
+    interfaces that do not display Appropriate Legal Notices, your
+    work need not make them do so.
+
+  A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit.  Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+  6. Conveying Non-Source Forms.
+
+  You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+    a) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by the
+    Corresponding Source fixed on a durable physical medium
+    customarily used for software interchange.
+
+    b) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by a
+    written offer, valid for at least three years and valid for as
+    long as you offer spare parts or customer support for that product
+    model, to give anyone who possesses the object code either (1) a
+    copy of the Corresponding Source for all the software in the
+    product that is covered by this License, on a durable physical
+    medium customarily used for software interchange, for a price no
+    more than your reasonable cost of physically performing this
+    conveying of source, or (2) access to copy the
+    Corresponding Source from a network server at no charge.
+
+    c) Convey individual copies of the object code with a copy of the
+    written offer to provide the Corresponding Source.  This
+    alternative is allowed only occasionally and noncommercially, and
+    only if you received the object code with such an offer, in accord
+    with subsection 6b.
+
+    d) Convey the object code by offering access from a designated
+    place (gratis or for a charge), and offer equivalent access to the
+    Corresponding Source in the same way through the same place at no
+    further charge.  You need not require recipients to copy the
+    Corresponding Source along with the object code.  If the place to
+    copy the object code is a network server, the Corresponding Source
+    may be on a different server (operated by you or a third party)
+    that supports equivalent copying facilities, provided you maintain
+    clear directions next to the object code saying where to find the
+    Corresponding Source.  Regardless of what server hosts the
+    Corresponding Source, you remain obligated to ensure that it is
+    available for as long as needed to satisfy these requirements.
+
+    e) Convey the object code using peer-to-peer transmission, provided
+    you inform other peers where the object code and Corresponding
+    Source of the work are being offered to the general public at no
+    charge under subsection 6d.
+
+  A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+  A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling.  In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage.  For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product.  A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+  "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source.  The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+  If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information.  But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+  The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed.  Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+  Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+  7. Additional Terms.
+
+  "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law.  If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+  When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it.  (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.)  You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+  Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+    a) Disclaiming warranty or limiting liability differently from the
+    terms of sections 15 and 16 of this License; or
+
+    b) Requiring preservation of specified reasonable legal notices or
+    author attributions in that material or in the Appropriate Legal
+    Notices displayed by works containing it; or
+
+    c) Prohibiting misrepresentation of the origin of that material, or
+    requiring that modified versions of such material be marked in
+    reasonable ways as different from the original version; or
+
+    d) Limiting the use for publicity purposes of names of licensors or
+    authors of the material; or
+
+    e) Declining to grant rights under trademark law for use of some
+    trade names, trademarks, or service marks; or
+
+    f) Requiring indemnification of licensors and authors of that
+    material by anyone who conveys the material (or modified versions of
+    it) with contractual assumptions of liability to the recipient, for
+    any liability that these contractual assumptions directly impose on
+    those licensors and authors.
+
+  All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10.  If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term.  If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+  If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+  Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+  8. Termination.
+
+  You may not propagate or modify a covered work except as expressly
+provided under this License.  Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+  However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+  Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+  Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License.  If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+  9. Acceptance Not Required for Having Copies.
+
+  You are not required to accept this License in order to receive or
+run a copy of the Program.  Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance.  However,
+nothing other than this License grants you permission to propagate or
+modify any covered work.  These actions infringe copyright if you do
+not accept this License.  Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+  10. Automatic Licensing of Downstream Recipients.
+
+  Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License.  You are not responsible
+for enforcing compliance by third parties with this License.
+
+  An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations.  If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+  You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License.  For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+  11. Patents.
+
+  A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based.  The
+work thus licensed is called the contributor's "contributor version".
+
+  A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version.  For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+  Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+  In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement).  To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+  If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients.  "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+  If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+  A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License.  You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+  Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+  12. No Surrender of Others' Freedom.
+
+  If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all.  For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+  13. Use with the GNU Affero General Public License.
+
+  Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work.  The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+  14. Revised Versions of this License.
+
+  The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+  Each version is given a distinguishing version number.  If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation.  If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+  If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+  Later license versions may give you additional or different
+permissions.  However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+  15. Disclaimer of Warranty.
+
+  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+  16. Limitation of Liability.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+  17. Interpretation of Sections 15 and 16.
+
+  If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+  If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+    <program>  Copyright (C) <year>  <name of author>
+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+  You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<https://www.gnu.org/licenses/>.
+
+  The GNU General Public License does not permit incorporating your program
+into proprietary programs.  If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.  But first, please read
+<https://www.gnu.org/licenses/why-not-lgpl.html>.

Armory/Armory/Makefile

@@ -0,0 +1,32 @@
+########################################################################
+#
+#
+#########################################################################
+
+CC	=	g++
+
+#COPTS	=	-fPIC -DLINUX -O2 -std=c++17 -lpthread
+COPTS	=	-fPIC -DLINUX -g -O0 -Wall -std=c++17 -lpthread
+
+ROOTLIBS = `root-config --cflags --glibs`
+
+ALL = Mapper EventBuilder#AnasenMS
+
+#########################################################################
+
+all	:	$(ALL)
+
+clean	:
+		/bin/rm -f $(OBJS) $(ALL)
+
+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)
+
+EventBuilder : EventBuilder.cpp ../ClassData.h fsuReader.h ../Hit.h
+		@echo "--------- making EventBuilder"
+		$(CC) $(COPTS) -o EventBuilder EventBuilder.cpp $(ROOTLIBS) 

Armory/Armory/Mapper.cpp

@@ -0,0 +1,182 @@
+#include <string>
+#include <cstdio>
+
+#include <TROOT.h>
+#include <TTree.h>
+#include <TFile.h>
+#include <TMath.h>
+#include <TBenchmark.h>
+
+#include "../mapping.h"
+#include "ClassDet.h"
+
+//===============================
+int main(int argc, char **argv){
+
+  printf("=========================================\n");
+  printf("===               Mapper              ===\n");
+  printf("=========================================\n");  
+  if (argc != 2)    {
+    printf("Incorrect number of arguments:\n");
+    printf("%s [inFile]\n", argv[0]);
+    printf("\n\n");
+    return 1;
+  }
+
+  ///============= read input
+  std::string inFileName = argv[1];
+
+  PrintMapping();
+
+  TFile * inFile = new TFile(inFileName.c_str(), "READ");
+  TTree * tree = (TTree*) inFile->Get("tree");
+  unsigned long long totnumEntry = tree->GetEntries();
+
+  ULong64_t       evID;
+  UInt_t          multi;
+  UShort_t         sn[MAXMULTI];
+  UShort_t         ch[MAXMULTI];
+  UShort_t          e[MAXMULTI];
+  UShort_t         e2[MAXMULTI];
+  ULong64_t       e_t[MAXMULTI];
+  UShort_t        e_f[MAXMULTI];
+
+  tree->SetBranchAddress("evID",  &evID);
+  tree->SetBranchAddress("multi", &multi);
+  tree->SetBranchAddress("sn",    sn);
+  tree->SetBranchAddress("ch",    ch);
+  tree->SetBranchAddress("e",     e);
+  tree->SetBranchAddress("e2",    e2);
+  tree->SetBranchAddress("e_t",   e_t);
+  tree->SetBranchAddress("e_f",   e_f);
+
+  ///================== new tree
+  TString outFileName = inFileName;
+  TString runStr = outFileName; 
+  int pos = outFileName.Last('/');
+  pos = outFileName.Index("_", pos+1); // find next "_"
+  runStr.Remove(0, pos+1);
+  runStr.Remove(3);
+  pos = outFileName.Index("_", pos+1); // find next "_"
+  outFileName.Remove(pos); // remove the rest
+  outFileName += "_mapped.root";
+
+  ULong_t eventID;
+  UInt_t run = runStr.Atoi();
+  
+  Det sx3;
+  Det qqq;
+  Det pc ;
+
+  printf(" Raw root file : %s\n", inFileName.c_str());
+  printf("           Run : %03d\n", run);
+  printf("   total Entry : %lld \n", totnumEntry);
+  printf(" Out file name : %s \n", outFileName.Data());
+
+  TFile * saveFile = new TFile( outFileName,"RECREATE");
+  TTree * newTree =  new TTree("tree","tree");
+  
+
+  newTree->Branch("evID",   &eventID,"eventID/l"); 
+  newTree->Branch("run",        &run,"run/i");  
+
+  newTree->Branch("sx3Multi", &sx3.multi, "sx3Multi/s");
+  newTree->Branch("sx3ID",    &sx3.id,    "sx3ID[sx3Multi]/s");
+  newTree->Branch("sx3Ch",    &sx3.ch,    "sx3Ch[sx3Multi]/s");
+  newTree->Branch("sx3E",     &sx3.e,     "sx3Energy[sx3Multi]/s");
+  newTree->Branch("sx3T",     &sx3.t,     "sx3Time[sx3Multi]/l");
+
+  newTree->Branch("qqqMulti", &qqq.multi, "qqqMulti/s");
+  newTree->Branch("qqqID",    &qqq.id,    "qqqID[qqqMulti]/s");
+  newTree->Branch("qqqCh",    &qqq.ch,    "qqqCh[qqqMulti]/s");
+  newTree->Branch("qqqE",     &qqq.e,     "qqqEnergy[qqqMulti]/s");
+  newTree->Branch("qqqT",     &qqq.t,     "qqqTime[qqqMulti]/l");
+  newTree->Branch("qqqSN",    &qqq.sn,     "qqqSN[qqqMulti]/s");
+
+  newTree->Branch("pcMulti", &pc.multi, "pcMulti/s");
+  newTree->Branch("pcID",    &pc.id,    "pcID[pcMulti]/s");
+  newTree->Branch("pcCh",    &pc.ch,    "pcCh[pcMulti]/s");
+  newTree->Branch("pcE",     &pc.e,     "pcEnergy[pcMulti]/s");
+  newTree->Branch("pcT",     &pc.t,     "pcTime[pcMulti]/l");
+
+  ///================== looping old tree and apply mapping
+
+  //clock   
+  // TBenchmark clock;
+  // Bool_t shown;
+
+  for( unsigned long long ev = 0; ev < totnumEntry; ev++){
+    tree->GetEntry(ev);
+
+    eventID = evID;
+    sx3.multi = 0;
+    qqq.multi = 0;
+    pc.multi = 0;
+
+    qqq.Clear();
+
+    for( unsigned int i = 0; i < multi; i++){
+
+      // printf("%10u/%10u| %5d, %2u, %6u, %14llu\n", i, multi, sn[i], ch[i], e[i], e_t[i] );
+      
+      //globalCh = digi-ID * nCh(digi-iD) + ch
+      int globalCh = -1;
+
+      for( int j = 0; j < nBd; j++){
+        if( board.at(j) == sn[i]){
+          globalCh = (sn[i] > 1000 ? j * 64 :  7*64 + (j-7) * 16) + ch[i]; //& = number V1740
+          break;
+        }
+      }
+      
+      if( globalCh == -1) printf("ev %llu\n", ev);
+
+      unsigned short ID = mapping[globalCh];
+
+      //=================================== sx3
+      if( ID < 10000 ) {
+        sx3.id[sx3.multi] = ID / 100;
+        sx3.ch[sx3.multi] = ID % 100;
+        sx3.e[sx3.multi] = e[i];
+        sx3.t[sx3.multi] = e_t[i];
+        sx3.multi ++;
+      }
+
+      //=================================== qqq
+      if( 10000 <= ID && ID < 20000 ) {
+        qqq.id[qqq.multi] = (ID - 10000) / 100;
+        qqq.ch[qqq.multi] = (ID - 10000) % 100;
+        qqq.e[qqq.multi] = e[i];
+        qqq.t[qqq.multi] = e_t[i];
+        qqq.sn[qqq.multi] = sn[i];
+        qqq.multi ++;
+      }
+
+      //=================================== pc
+      if( 20000 <= ID && ID < 30000 ) {
+        pc.id[pc.multi] = (ID - 20000) / 100;
+        pc.ch[pc.multi] = (ID - 20000) % 100;
+        pc.e[pc.multi] = e[i];
+        pc.t[pc.multi] = e_t[i];
+        pc.multi ++;
+      }
+    }
+    
+    saveFile->cd(); //set focus on this file
+    newTree->Fill();
+    
+    if( eventID % 100 == 0 ) printf("%6lu/%6llu [%2d%%]\n\033[A\r", eventID, totnumEntry, TMath::Nint((eventID+1)*100./totnumEntry));
+
+  }
+
+  inFile->Close();
+
+  saveFile->cd(); //set focus on this file
+  newTree->Write(); 
+  UInt_t eventNumber = newTree->GetEntries();
+  saveFile->Close();
+  printf("-------------- done,  %u\n", eventNumber);
+  
+  return 0;
+
+}

Armory/Armory/anasenMS.cpp

@@ -0,0 +1,280 @@
+#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"
+
+//======== Gerneate light particle based on reaction
+// find out the CalTrack and the real track
+// find out the Q-value uncertaintly
+
+int main(int argc, char **argv){
+
+  printf("=========================================\n");
+  printf("===       ANASEN Monte Carlo          ===\n");
+  printf("=========================================\n");  
+  
+  int numEvent = 1000000;
+  if( argc >= 2 ) numEvent = atoi(argv[1]);
+
+  //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",numEvent);
+  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;
+
+  //================================= Calculate event
+  for( int i = 0; i < numEvent ; 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();
+
+    }
+
+    tree->Fill();
+
+    //#################################################################### Timer  
+    clock.Stop("timer");
+    Double_t time = clock.GetRealTime("timer");
+    clock.Start("timer");
+
+    if ( !shown ) {
+      if (fmod(time, 10) < 1 ){
+        printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60);
+        shown = 1;
+      }
+    }else{
+      if (fmod(time, 10) > 9 ){
+        shown = 0;
+      }
+    }
+
+  }
+
+  tree->Write();
+  int count = tree->GetEntries();
+  saveFile->Close();
+
+  printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
+
+  delete anasen;
+
+  return 0;
+
+}

Armory/Armory/constant.h

@@ -0,0 +1,105 @@
+/***********************************************************************
+ * 
+ *  This is constant.h, to provide various physical constants.
+ * 
+ *-------------------------------------------------------
+ *  created by Ryan (Tsz Leung) Tang, Nov-18, 2018
+ *  email: goluckyryan@gmail.com
+ * ********************************************************************/
+
+#ifndef constant
+#define constant
+#include <cmath>
+
+const double pi   = acos(-1.0);
+const double E    = 2.718281828459 ;
+const double hbar_SI = 1.054571628e-34;   //Js
+const double kB   = 1.3806504e-23;     //JK^-1
+const double e    = 1.602176487e-19;       //C
+const double c_SI    = 299792458;         //ms^-1
+const double me_SI   = 9.10938215e-31 ;   //kg
+const double mp_SI   = 1.672621637e-27 ;  //kg
+const double mn_SI   = 1.67492729e-27 ;   //kg
+const double NA   = 6.022141e+23 ;     //mol^-1
+
+const double deg2rad = pi/180 ;
+const double rad2deg = 180/pi ;
+
+//======================================================================
+const double amu  = 931.49432; // MeV/c^2
+const double hbarc = 197.326979; // MeV fm;
+const double c = 299.792458; // mm/ns; 
+const double ee = 1.439964454; // MeV.fm
+
+//======================================================================
+double kg2MeV(double m){
+    return m*c_SI*c_SI/e/1e6;
+}
+
+double T2Brho(double mass, int Z, int A, double T){
+    //mass in MeV
+    // Z in e
+    // T in MeV/A
+    double gamma = (T*A + mass)/mass;
+    double beta = sqrt(1-1/gamma/gamma);
+    return mass*beta*gamma/Z/c;
+}
+
+double Brho2T(double mass, int Z, int A, double Brho){
+    //mass in MeV
+    // Z in e
+    return (sqrt(pow(Brho*Z*c,2)+mass*mass)-mass)/A;
+}
+
+double T2beta(double mass, int A, double T){
+	double gamma = 1.0 + T*A/mass;
+	return sqrt(1-1/gamma/gamma);
+}
+
+double ev2nm(double eV){
+  // photon energy to nm
+  return hbarc/2/pi/eV;
+}
+
+//======================================================================
+const double mp = kg2MeV(mp_SI);
+const double mn = kg2MeV(mn_SI);
+const double hbar = 197.326979;
+
+
+//======================================================================
+inline std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){
+
+  std::vector<std::string> output;
+
+  size_t pos;
+  do{
+    pos = tempLine.find(splitter); /// fine splitter
+    if( pos == 0 ){ ///check if it is splitter again
+      tempLine = tempLine.substr(pos+1);
+      continue;
+    }
+
+    std::string secStr;
+    if( pos == std::string::npos ){
+      secStr = tempLine;
+    }else{
+      secStr = tempLine.substr(0, pos+shift);
+      tempLine = tempLine.substr(pos+shift);
+    }
+
+    ///check if secStr is begin with space
+    while( secStr.substr(0, 1) == " ") secStr = secStr.substr(1);
+    
+    ///check if secStr is end with space
+    while( secStr.back() == ' ') secStr = secStr.substr(0, secStr.size()-1);
+
+    output.push_back(secStr);
+    ///printf(" |%s---\n", secStr.c_str());
+    
+  }while(pos != std::string::npos );
+
+  return output;
+}
+
+#endif

Armory/Makefile

@@ -27,6 +27,6 @@ Mapper	:	Mapper.cpp  ../mapping.h ClassDet.h
 # 		@echo "--------- making ANASEN Monte Carlo"
 # 		$(CC) $(COPTS) -o AnasenMS anasenMS.cpp $(ROOTLIBS)
 
-EventBuilder : EventBuilder.cpp ClassData.h fsuReader.h Hit.h
+EventBuilder : EventBuilder.cpp ../ClassData.h fsuReader.h ../Hit.h
 		@echo "--------- making EventBuilder"
 		$(CC) $(COPTS) -o EventBuilder EventBuilder.cpp $(ROOTLIBS) 

MakeVertex.C

@@ -39,6 +39,7 @@ bool realtime = true;
 const double source_vertex = 53; //53
 const double qqq_z = 100.0;
 const double anode_gain = 1.5146e-5; //channels --> MeV
+std::string dataset = "26Al"; //"17F"
 
 TApplication *app=NULL;
 TH1F *hha=NULL,*hhc=NULL;
@@ -124,7 +125,7 @@ void MakeVertex::Begin(TTree * /*tree*/)
     }
 
     // Load PC Calibrations
-    std::ifstream inputFile("slope_intercept_results.txt");
+    std::ifstream inputFile("slope_intercept_results_"+dataset+".txt");
     if (inputFile.is_open())
     {
         std::string line;
@@ -183,7 +184,7 @@ void MakeVertex::Begin(TTree * /*tree*/)
     }
 
     {
-        std::ifstream infile("sx3cal/17F/backgains.dat");
+        std::ifstream infile("sx3cal/"+dataset+"/backgains.dat");
         std::string temp;
         int backpos, frontpos, clkpos;
         if (infile.is_open())
@@ -191,13 +192,13 @@ void MakeVertex::Begin(TTree * /*tree*/)
                 ;//std::cout << sx3BackGain[clkpos][frontpos][backpos] << std::endl;
         infile.close();
 
-        infile.open("sx3cal/17F/frontgains.dat");
+        infile.open("sx3cal/"+dataset+"/frontgains.dat");
         if (infile.is_open())
             while(infile>>clkpos>>temp>>temp>>frontpos>>sx3FrontOffset[clkpos][frontpos]>>sx3FrontGain[clkpos][frontpos])
                 ;//std::cout << sx3FrontOffset[clkpos][frontpos] << " " << sx3FrontGain[clkpos][frontpos] << std::endl;
         infile.close();
 
-        infile.open("sx3cal/17F/rightgains.dat");
+        infile.open("sx3cal/"+dataset+"/rightgains.dat");
         if (infile.is_open())
             while(infile>>clkpos>>frontpos>>temp>>sx3RightGain[clkpos][frontpos]) {
                 sx3RightGain[clkpos][frontpos]=TMath::Abs(sx3RightGain[clkpos][frontpos]);
@@ -221,7 +222,6 @@ void MakeVertex::Begin(TTree * /*tree*/)
         can1->Modified();
         can1->Update();
         can1->BuildLegend();
-
         can2->cd();
         frame->Draw();
         for(int i=0; i<24; i++) {
@@ -291,11 +291,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
     qqq.CalIndex();
     pc.CalIndex();
 
-/*    for (int i = 0; i < pc.multi; i++)
+    std::vector<Event> sx3Events;
-    {
-        std::cout << pc.index[i] << " " << pc.e[i] << " " << std::endl;
-    }
-*/    std::vector<Event> sx3Events;
     if(sx3.multi>1) {
         std::array<sx3det,24> Fsx3;
         //std::cout << "-----" << std::endl;
@@ -364,21 +360,6 @@ Bool_t MakeVertex::Process(Long64_t entry)
     int qqqCount = 0;
     int qqqAdjCh = 0;
     // REMOVE WHEN RERUNNING USING THE NEW CALIBRATION FILE
-    for (int i = 0; i < qqq.multi; i++)
-    {
-        //if ((qqq.id[i] == 3 || qqq.id[i] == 1) && qqq.ch[i] < 16)
-        if (qqq.id[i] == 1 && qqq.ch[i] < 16) //for run 12, 26Al
-        {
-            qqq.ch[i] = 16 - qqq.ch[i];
-        }
-    }
-    for (int i = 0; i < qqq.multi; i++)
-    {
-        if (qqq.id[i] == 0 && qqq.ch[i] >= 16)
-        {
-        }
-    }
-
     std::vector<Event> QQQ_Events, PC_Events;
     std::vector<Event> QQQ_Events_Raw, PC_Events_Raw;
     std::vector<Event> QQQ_Events2; //clustering done
@@ -519,7 +500,8 @@ Bool_t MakeVertex::Process(Long64_t entry)
                     double x = rho * TMath::Cos(theta);
                     double y = rho * TMath::Sin(theta);
                     hitPos.SetXYZ(x, y, qqq_z);
-                    if(realtime) qqqg->SetPoint(0,hitPos.X(),hitPos.Y(),hitPos.Z());
+                    //if(realtime) qqqg->SetPoint(0,hitPos.X(),hitPos.Y(),hitPos.Z());
+                    if(realtime) qqqg->AddPoint(hitPos.X(),hitPos.Y(),hitPos.Z());
                     qqqenergy = eRingMeV;
                     qqqtimestamp = tRing;
                     HitNonZero = true;
@@ -572,6 +554,10 @@ Bool_t MakeVertex::Process(Long64_t entry)
         {
             cathodeT = static_cast<double>(pc.t[i]);
             cathodeIndex = pc.index[i] - 24;
+
+            //int flipped_index = 23-cathodeIndex;
+            //cWireEvents[flipped_index] = std::tuple(flipped_index,pc.e[i],static_cast<double>(pc.t[i]));
+
             cWireEvents[pc.index[i]-24] = std::tuple(pc.index[i]-24,pc.e[i],static_cast<double>(pc.t[i]));
             if(realtime) hhc->SetBinContent(hhc->FindFixBin(cathodeIndex),pc.e[i]);
         }
@@ -702,6 +688,7 @@ Bool_t MakeVertex::Process(Long64_t entry)
             if(PCSX3TimeCut) {
                 plotter->Fill2D("xyplot_sx3"+std::to_string(sx3event.ch2/4),100,-100,100,100,-100,100,sx3event.pos.X(),sx3event.pos.Y());
                 plotter->Fill2D("xyplot_sx3"+std::to_string(sx3event.ch2/4),100,-100,100,100,-100,100,pcevent.pos.X(),pcevent.pos.Y());
+                plotter->Fill2D("pcz_vs_pcphi_TimeCut",600,-200,200,120,-360,360,pcevent.pos.Z(),pcevent.pos.Phi()*180/M_PI); //x-axis is all Si det, y-axis is PC anode+cathode only
             }
             double sx3rho = 88.0;//approximate barrel radius
             double sx3z = sx3event.pos.Z()+(75.0/2.0)-3.0; //w.r.t target origin at 90 for run12
@@ -898,7 +885,8 @@ Bool_t MakeVertex::Process(Long64_t entry)
         else {
             anodeIntersection = TVector3(x, y, z);
             if(realtime) {
-                crossoverg->SetPoint(0,x,y,z);
+                //crossoverg->SetPoint(0,x,y,z);
+                crossoverg->AddPoint(x,y,z);
             }
             //std::cout << "Anode Intersection: " << anodeIntersection.X() << ", " << anodeIntersection.Y() << ", " << anodeIntersection.Z() << " " << aIDMax <<  std::endl;
         }
@@ -910,8 +898,8 @@ Bool_t MakeVertex::Process(Long64_t entry)
     }
 
     if(anodeIndex!=-1 && cathodeIndex !=-1 && hitPos.Perp()!=0 && anodeIntersection.Perp()!=0 && realtime && PCQQQPhiCut && PCQQQTimeCut) {
-        can1->Modified();
+        //can1->Modified();
-        can1->Update();
+        //can1->Update();
         TVector3 x2(anodeIntersection);
         TVector3 x1(hitPos);
         TVector3 v = x2-x1;
@@ -930,8 +918,9 @@ Bool_t MakeVertex::Process(Long64_t entry)
             //pla[anodeW.first]->SetLineStyle(kLine);
         }
         //can2->Modified();
-        can2->Update();
+        //can2->Update();
-        while(can1->WaitPrimitive());
+        //while(can1->WaitPrimitive());
+
         //pla[anodeIndex]->SetLineWidth(1);
         //pla[anodeIndex]->SetLineStyle(kDotted);
         for(auto anodeW: anodeHits) {
@@ -942,7 +931,6 @@ Bool_t MakeVertex::Process(Long64_t entry)
             plc[cathodeIndex]->SetLineStyle(kDotted);
             plc[cath.first]->SetLineWidth(1);
         }
-
     }
 
     //for (double Tz = 60; Tz <= 100; Tz += 1.0)
@@ -1226,4 +1214,10 @@ Bool_t MakeVertex::Process(Long64_t entry)
 void MakeVertex::Terminate()
 {
     plotter->FlushToDisk();
+    can1->Modified();
+    can1->Update();
+    can2->Modified();
+    can2->Update();
+    while(can1->WaitPrimitive());
+    while(can2->WaitPrimitive());
 }

eloss_calculations/make_eloss_table.C

@@ -0,0 +1,13 @@
+#include "/home/sud/Desktop/Software2/propagator/elastcaller.h"
+void make_eloss_table() {
+    double einput = 6.0, estepnow;
+    double target_thickness_unit = 1e-3; //mg/cm2. 
+    double density = 0.0711;//mg/cm3
+    long i=0;
+    while(einput > 4.0) {
+        std::cout << "After " << i << " steps, 4He is at " << einput << " MeV after penetrating " << i*target_thickness_unit << " mg/cm2 " << i*target_thickness_unit/density << " cm of HeCO2" << std::endl;
+        estepnow = slowmedown("4He",einput,"3(12C)6(16O)97(4He)",target_thickness_unit);
+        einput = estepnow;
+        i+=1;
+    }
+}

eloss_calculations/stopit/desorb.for

@@ -0,0 +1,811 @@
+	subroutine desorb(ianz,zp,ap,ep,loste)
+
+C **  CALCULATES ENERGY LOSS IN AN ABSORBER SANDWICH *******
+C ****  ENERGY DEPOSIT IN SECTIONS OF IONIZATION  **********
+C ***************  CHAMBER (DE1, DE2, AND DE3)  ************
+
+
+	PARAMETER (DSEP0=4.0,isold0=1)
+
+        COMMON ISG(19),INN(19),DEN(19),THK(19),PRS(19),XLN(19),ARDEN(19)
+	1      ,ZNUMB(19,4),ANUMB(19,4),ELNUM(19,4),CONCN(19,4)
+	1         ,PRDEN(19,4),PRTHK(19,4)
+	DIMENSION ZNUMBW(4),ANUMBW(4),ELNUMW(4),CONCNW(4)
+	1         ,PRDENW(4),PRTHKW(4)
+	DIMENSION E(19),DE(19),xmem(19)
+	DIMENSION TOUT(19),TOUTE(19)
+	dimension eptable(50,500,2),emintabz(50),ianzv(5)
+	real loste(19)
+	DATA io1,IO2,IO3,io0/9,11,12,1/
+	data iopt,ianzi,ianzide/1,2,1/
+c
+	save eltimo,izmax,izth
+
+c       the mass table is to be used only for iopt = 5,6
+c       use atomic masses to average for isotipic composition.
+c       taken from Formulas Facts and Constants, H. J. Fischbeck and
+c       K. H. Fischbeck. Springer - Verlag 1987 2nd ed, pages 164-183.
+
+        dimension amass(70)
+        data amass/1.01,4.00,6.94,9.01,10.81,12.01,14.01,16.00,19.00
+     1,20.18,22.99,24.31,26.98,28.09,30.97,32.07,35.45,39.95
+     2,39.10,40.08,44.96,47.88,50.94,52.00,54.94,55.85,58.93
+     3,58.69,63.55,65.39,69.72,72.59,74.92,78.96,79.90,83.80
+     4,85.47,87.62,88.91,91.22,92.91,95.94,98.,101.07,102.91
+     5,106.42,107.87,112.41,114.82,118.71,121.75,127.60,126.90
+     6,131.29,132.91,137.33,138.91,140.12,140.91,144.24,147.
+     7,150.36,151.96,157.25,158.93,162.5,164.93,167.26,168.93
+     8,173.04/
+
+c       for Z > 70, you are in trouble !!
+c
+c      open(unit=IO1, status='OLD', file='absorb.inp')
+      open(unit=IO2, status='OLD', file='desorb.out')
+c      rewind IO1
+c      rewind IO2
+c
+10	CONTINUE
+
+
+C      IOPT = 1 - SUPPLY ENERGY OF PARTICLE ENTERING
+C                 THE ABSORBER ARRAY AND GET LOSS AND
+C                 RANGES
+C      IOPT = 2 - SUPPLY TARGET, PROECTILE AND EJECTILE
+C                 INFO. AND THEN GO THROUGH ABSORBER
+C                 SANDWICH
+C      IOPT = 3 - CALCULATE ENERGY DEPOSITS IN DETECTOR
+C                 DETECTOR DIMENSIONS ARE STANDARD AND
+C                 THE VARIABLE -'IDET' - CHOOSES BETWEEN
+C                 VARIETY OF AVAILABLE DETECTORS
+C      IOPT = 4 - FINDS MAXIMUM ENERGY THAT CAN BE STOPPED IN
+C                 IANZ ELEMENTS OF THE SANDWICH FOR GIVEN
+C                 ZP, AP.
+C                 WHEN CALCULATION IS FINISHED, THE PROGRAM READS
+C                 IN NEW VALUES OF ZP, AP AND RESTARTS. TO END
+C                 THE PROGRAM, GIVE ZP < 0.
+C                 IN ORDER TO HELP THE SPEED OF THE PROGRAM,
+C                 GIVE THE PARTICLE'S "Z" IN increasing  ORDER.
+C      IOPT = 5 - STORES ARRAYS OF Edet AS A FUNCTION OF INCIDENT    
+C                 ENERGY AND THE PARTICLE'S ID (Z,A)
+C                 ARRAY LABELED  eptable(Z-Zth,Einc,ipunch)
+C                 ipunch = 1 stopped,  = 2 punched through
+c                 Einc = E(incident)/detable
+C                 Zth  = lowest Z considered - 1
+C
+C ************************************************************
+
+c	read(io1,*) iopt
+c1	FORMAT(10I)
+	DSEP = DSEP0
+	isold = isold0
+c
+	IF(iopt.LT.0) GO TO 2000
+
+
+C ***********************************************************
+c
+	if(iopt.eq.3) then
+      open(unit=IO3, status='OLD', file='absorbv.out')
+      rewind IO3
+	endif
+c
+	if(iopt.ge.5) then
+		do itpun = 1,2
+		do itblz = 1,50
+		emintabz(itblz) = 0.
+		do itble = 1,500
+		eptable ( itblz,itble,itpun ) = 0.
+		enddo
+		enddo
+		enddo
+c
+		open(unit=io0, status='UNKNOWN', file='abs.tbl')
+		rewind io0
+c
+	else
+	endif
+
+C ***********************************************************
+
+c      ianz  = number of elements in absorber "sandwich" - the 
+c              particle deposits all its energy in these layers.
+c      ianzi = index of last layer in which the energy of the 
+c              particle is not recorded - this unreecorded energy
+c              is used in the DST production in two modes:
+c          when making DST tape from data:
+c              Since the detector records only deposited energy
+c              the output tables are used to correct for this 
+c              deficinecy and for a given dharge and mass extrapolate
+c              the measured energy to the target exit energy
+c          when making DST tape from model calculations:
+c              The "lost" energy is a source of broadening of the 
+c              energy spectra due to straggling - this "smearing" 
+c              is estimated and superposed on the calculated spectra.
+c	ianzide = element # for DE calculation
+c
+c	read(IO1,*) ianz,ianzi,ianzide
+c
+c	ianzv(5)
+c	     = Index of layer where exiting particle velocit is 
+c	       calculated (only for option 3) max = 5
+c
+c	if(iopt.eq.3) read (IO1,*) knz
+c	if(iopt.eq.3) read (IO1,*) (ianzv(k),k=1,knz)
+C ***********************************************************
+
+C      AP = PROJECTILE MASS
+C      ZP = PROJECTILE CHARGE
+C      EP = PROJECTILE ENERGY
+c
+c	read(io1,*) zp,ap,ep
+2	FORMAT(f10.3)
+
+C *************************************************************
+c
+c      zth= threshold Z for incident energy tble calc. (iopt=5,6)
+c      zmax = maximum z for table calculation
+c      detable = the energy step size used for array storage   
+c      emin = starting incident energy for table
+c      emax = mazximum incident energy for table calculation
+c	EP is ignored when iopt = 5 or 6
+c
+	if(iopt.ge.5) then
+		read(IO1,*) zth,zmax,emintab,emaxtab,detable
+		eltimeo = secnds(0.0)  !  start timing
+		izth = ifix (zth + 0.01)
+		izmax = ifix (zmax + 0.01)
+c
+c	if detable > emintab change emintab to 1/2 * detable
+c
+		if(emintab.lt.detable) emintab = 0.5*detable
+		write(io2,291) zth+1,zmax,emintab,emaxtab,detable
+	else
+	endif
+c
+C ************************************************************
+
+C      IN THE FOLLOWING THE LISTED VARIABLES ARE INDEXED
+C      THE INDICES I AND J STAND FOR THE FOLLOWING:
+C      - I -  SERIAL NUMBER OF ABSORBER LAYER (<20)
+C      - J -  SERIAL NUMBER OF ELEMENT WITHIN LAYER (<5)
+
+C *************************************************************
+
+
+C       ISG(I) = 0  - FOR SOLID ABSORBERS
+C          (I) = 1  - FOR GASEOUS ABSORBERS
+C       INN(I) = NUMBER OF ELEMENTS IN ONE LAYER
+C             (E.G. CH4 HAS TWO ELEMENTS C AND H)
+C       DEN(I) = DENSITY OF ABSORBE (FOR SOLIDS)
+C       THK(I) = THICKNESS OF ABSORBER IN MG/CMSQ (FOR SOLIDS)
+C       PRS(I) = PRESSURE (IN MM HG) FOR GAS ABSORBER
+C       XLN(I) = PHYSICAL LENGTH OF ABSORBER ( IN CM ) FOR GAS
+
+C  **** ABSORBER COMPOSITION ****
+
+C      ELNUM(I,J) = NUMBER OF ATOMS OF ELEMENT J IN LAYER I
+C      CONCN(I,J) = CONCENTRATION (MOLAR) OF ELEMENT J IN LAYER I
+C      ANUMB(I,J) = MASS   NUMBER   OF ELEMENT J IN LAYER I
+C      ZNUMB(I,J) = ATOMIC NUMBER   OF ELEMENT J IN LAYER I
+
+C      PRDEN(I,J) = PARTIAL DENSITY OF ELEMENT J IN LAYER I
+C      PRTHK(I,J) = PARTIAL THICKNESS OF ELEMENT J IN LAYER I (MG/CMSQ)
+
+C      MAXIMUM OF NINETEEN LAYERS SPECIFIED
+
+
+C *************************************************************
+
+	DO 100 ISN=1,IANZ
+c		read(io1,*) ISG(ISN),INN(ISN)
+		IF(ISG(ISN).EQ.1) GO TO 50
+c		read(io1,*) DEN(ISN),THK(ISN)
+		TOUT(ISN)=THK(ISN)/(DEN(ISN)*1000.)
+		TOUTE(ISN)=TOUT(ISN)/2.54
+		DO 20 IMN=1,INN(ISN)
+c			read(io1,*) ANUMB(ISN,IMN),ZNUMB(ISN,IMN),
+c	1		ELNUM(ISN,IMN),CONCN(ISN,IMN)
+20		CONTINUE
+		GO TO 100
+50              continue
+c		read(io1,*) PRS(ISN),XLN(ISN)
+		DO 60 IMN=1,INN(ISN)
+c			read(io1,2) ANUMB(ISN,IMN),ZNUMB(ISN,IMN),
+c	1		ELNUM(ISN,IMN),CONCN(ISN,IMN)
+60		CONTINUE
+100	CONTINUE
+
+C ****************************************************************
+
+	if(iopt.ne.5.and.iopt.ne.6) WRITE(IO2,101) ap,zp,ep
+	if(iopt.ne.5.and.iopt.ne.6) WRITE(IO2,102) ianz
+
+C *****************************************************************
+
+	DO 200 I=1,ianz
+		INNW=INN(I)
+		DO 210 J=1,INNW
+			ANUMBW(J)=ANUMB(I,J)
+			ZNUMBW(J)=ZNUMB(I,J)
+			ELNUMW(J)=ELNUM(I,J)
+			CONCNW(J)=CONCN(I,J)
+210		CONTINUE
+		DENW=DEN(I)
+		XLNW=XLN(I)
+		PRSW=PRS(I)
+		THKW=THK(I)
+		IF(ISG(I).EQ.1) GO TO 250
+		CALL SETABS(INNW,ANUMBW,ZNUMBW,ELNUMW,PRTHKW,THKW,PRDENW,DENW)
+		DO 230 J=1,INNW
+			PRDEN(I,J)=PRDENW(J)
+			PRTHK(I,J)=PRTHKW(J)
+230		CONTINUE
+		GO TO 200
+
+250		CALL SETABG(INNW,ANUMBW,ZNUMBW,ELNUMW,CONCNW,PRTHKW,THKW
+	1,	PRDENW,DENW,PRSW,XLNW)
+		DEN(I)=0.
+		THK(I)=0.
+		DO 270 J=1,INNW
+			PRDEN(I,J)=PRDENW(J)
+			PRTHK(I,J)=PRTHKW(J)
+			DEN(I)=DEN(I)+PRDEN(I,J)
+			THK(I)=THK(I)+PRTHK(I,J)
+270		CONTINUE
+200	CONTINUE
+
+C *************************************************************
+
+C      START CALCULATION AND DETAILED PRINTOUT
+
+C *************************************************************
+c
+	if(iopt.ge.5) then
+	ep = emintab
+	zp = zth + 1.
+	indexz = ifix (zp - zth + 0.001)
+	endif
+c
+299	continue   !  come here for new particle (zp change)
+
+	izp = ifix (zp+0.001)
+c
+	if(iopt.ge.5) then
+		if (izp.gt.70) then
+			write (6,*) 'no mass for Z = ',izp
+			stop
+		else
+			ap = amass(izp)
+
+c       The trick!. To calculate energy losses for deuterons and tritons,
+c       enter zth = 0.2 and 0.3 respectively.    E. Chavez jul/92
+
+			if (izp.eq.1) then
+				iap = ifix (zth*10.0 + 0.1)
+				ap = float (iap)
+			end if
+		end if
+	end if
+c
+	if (iopt.eq.6) then
+		ideltalay = 8                 ! choose DE3 for eloss signal
+		if(izp.eq.1) ideltalay = 16   ! choose DEh for eloss signal
+	endif
+c
+300	CONTINUE    ! come here for new energy (ep)
+c
+	EI=ep
+	XUPDN=-1.
+	EPS=0.0001
+	I1STPASS = 1
+
+	IF (iopt.EQ.4) GO TO 600
+
+	ipunch = 2
+	DO 510 I=1,IANZ    ! begin loop over absorber layers
+c
+	if(iopt.ge.5) go to 504
+	IF(ISG(I).EQ.0) WRITE(IO2,311) I,THK(I),TOUT(I),TOUTE(I),DEN(I)
+	IF(ISG(I).EQ.1) WRITE(IO2,312) I,THK(I),PRS(I),XLN(I),DEN(I)
+	DO 320 J=1,INN(I)
+	WRITE(IO2,321) ANUMB(I,J),ZNUMB(I,J),PRTHK(I,J)
+320	CONTINUE
+504	continue
+c
+c       XNS - initial no. of intervals for integration of DE
+		XNS = 2.
+c	EI = energy in   
+		CALL ADS(I,XUPDN,XNS,EPS,ap,zp,EI,DEI,ISTAT)
+c       DEI = energy out - energy in ( < 0. for energy loss)
+		EIOLD=EI
+c       E(I) = energy left after I'th element (EP-DE(1)-DE(2)+...)
+c       if particle stopped in detector this is equal to energy lost
+c       in remaining layers
+		DE(I) = DEI
+		E(I)  = EI + DEI
+		EI    = E(I)
+		INS=IFIX(XNS+0.001)
+c
+	if (iopt.ge.5) go to 505
+	WRITE(IO2,401) INS,EIOLD,EI
+	loste(i)=-1.*de(i)
+	if (EI.LT.EPS.OR.ISTAT.EQ.-1) WRITE(IO2,402) I
+505	continue
+        loste(ianz + 1)=e(ianz)
+c
+c	if particle stopped in layer beyond ianzi we must 
+c	check if iopt=5 or6 and calculate the energy loss in the
+c	front part (layers 1 thru ianzi).
+c
+	istore = I
+	if (EI.lt.EPS) ipunch=1
+c	control loop exit
+c	exit when particle runs out of energy in last layer
+	if (EI.LT.EPS.OR.ISTAT.EQ.-1) go to 701
+c	if interested in DE signal only (iopt=6) exit after 
+c	layer for which DE is sought was traversed
+	if(iopt.eq.6.and.I.ge.ianzide) go to 701
+c
+510	CONTINUE     ! end loop over absorber layers
+
+c	this part for iopt=5,6 stores incident energy values
+
+701	continue
+	if(iopt.ne.5.and.iopt.ne.6) go to 520
+c
+c  establish higher energy cutoof for next step with higher z
+c  should save time in calulating energies of particles stopped
+c  in the dead layer
+	if( I.le.ianzi) emintabz(izp)=ep
+c
+c  evalue = energy of particle when entering sensitive volume of det.
+c           when particle is stopped (ipunch=1) this is what is left
+c           once you take off the energy lost in the dead layer.
+	evalue = E(ianzi) 
+c         = when particle punches thouough detector its energy at the
+c           end is EI - subtract this from what it entered with (after
+c           dead layers) and again you got the energy deposited.
+	if(ipunch.eq.2) evalue = E(ianzi) - EI 
+c      roundoff errors could resutl in negative energies !!
+	if (evalue.lt.0.0) evalue = 0.0
+
+c  EI is current energy after last layer - is nonzero when particle
+c  punched through and to get energy deposited must subtract this 
+c  "left over" energy from the energy of particle had when it entered
+c  the detector's sensitive volume
+	indexe = ifix((ep + 0.001)/detable) + 1
+        indexz = ifix(zp - zth + 0.001)
+	if(iopt.eq.5) eptable(indexz,indexe,ipunch) = evalue
+	if(iopt.eq.6) ipunch = 1
+	if(iopt.eq.6) eptable(indexz,indexe,ipunch) = -DE(ianzide)
+
+c	now repeat calculation for same Z but new energy
+	ep = ep + detable 
+	if(ep.gt.emaxtab) go to 709
+	go to 300
+
+709	continue
+c	reset energy to emintab and up zp by one until we top
+c	zmax - this portion controls looping over Z!
+
+	do ieps = 1,ianz
+		E(ieps) = 0.
+		DE(ieps) = 0.
+	enddo
+
+	ep = emintabz(izp)
+	zp = zp + 1.
+	izp = ifix(zp + 0.001)
+	emintabz(izp) = emintabz(izp-1)
+
+	eltime = secnds (eltimeo)
+	itminutes = ifix(eltime/60.)
+	tminutes = float(itminutes)
+	tseconds = eltime - tminutes*60.
+	eltimeo = eltimeo + eltime
+	zpm1 = zp - 1.
+	if((izp-1).le.izmax) write(io2,714) zpm1,ap,tminutes,tseconds
+  714 format(' finished Z=',f3.0,'   A=',f4.0,' -  ',f5.0,
+     +' minutes and ',f3.0,' seconds elapsed')
+
+	if (izp.ge.izmax) go to 711
+	go to 299
+
+711	continue
+
+c	get here when iopt = 5 or 6 calculation is done
+c	now ready to store array on disk
+	write(io0,712) zth,zmax,emintab,emaxtab,detable
+
+	iemintab = ifix( (emintab + 0.001)/detable ) + 1
+	iemaxtab = ifix( (emaxtab + 0.001)/detable ) 
+	iztop = ifix(zp - 1. - zth + 0.001)
+
+	do ipp = 1,2  
+	do indexz = 1, iztop
+	iztab = indexz + ifix(zth+0.01)
+	imasstab = ifix(amass(iztab) + 0.1)
+	if(iztop.eq.1) imasstab = ifix(ap + 0.1)
+	write(io0,7712) indexz,iztab,imasstab,iemintab,iemaxtab
+7712	format(5i20)
+	itblow = iemintab 
+	do indexe = iemintab,iemaxtab,10
+	itbup = itblow + 9
+	write(io0,713) (eptable(indexz,itbe,ipp),itbe=itblow,itbup)
+	itblow = itbup + 1
+	enddo
+	enddo
+	enddo
+712	format(5e16.8)
+713	format(10e16.8)
+c
+	close (unit = io0)
+c
+C      MORE INPUT FOR NEW CALCULATION WITH SAME ABSORBERS
+
+520	CONTINUE
+
+
+c	read(io1,*) zp,ap,ep
+        zp = -1.
+	IF(zp.le.0.) GO TO 2000
+	izp = ifix(zp + 0.001)
+	if(ap.le.0.) AP = amass(izp)
+	IF(zp.gt.0.) WRITE(IO2,101) ap,zp,ep
+	GO TO 299
+
+600	DO I = 1, ianz
+		ILAY = I
+		XNS = 2.0
+		CALL ADS(I,XUPDN,XNS,EPS,ap,zp,EI,DEI,ISTAT)
+		EIOLD = EI
+		DE(I) = DEI
+		E(I) = EI + DEI
+c       E(I) = energy left after I'th element (EP+DE(1)+DE(2)+...)
+c       if particle stopped in detector this is equal to energy lost
+c       in remaining layers
+		xmem(i) = E(I)
+		EI = E(I)
+		INS = IFIX(XNS + 0.001)
+c		WRITE (IO2,613) ILAY,INS,EI,ISTAT
+		IF (EI.LE.0.0.OR.ISTAT.EQ.-1) GO TO 601
+	END DO
+
+
+601	IF (ISTAT.EQ.0)THEN
+		IF (EI.LT.0.003.AND.EI.GE.0.0) THEN
+	WRITE(IO2,611) ap,zp,ep,ILAY,xmem(5),xmem(6),xmem(7)
+		read(io1,*) zp,ap
+		izp = ifix (zp + 0.001)
+		if(ap.le.0.) ap = amass(izp)
+			IF (zp.LT.0.0) GO TO 2000
+			IPASS = 0
+			I1STPASS = 1
+			EI = ep
+			GO TO 600
+		END IF
+		isign = isold
+		isold = 1
+	ELSE
+		isign = - isold
+		isold = -1
+	END IF
+
+	IF (I1STPASS.gt.0) THEN
+		isign = 1
+		I1STPASS = 0
+		IF (ISTAT.EQ.0) THEN
+			DSEP =  - DSEP0
+		ELSE
+			DSEP = DSEP0
+		END IF
+	END IF
+
+C	IF THE INITIAL ENERGY WAS TOO LARGE, THEN THE ION WILL PUNCH THROUGH
+C	THE DETECTOR A NUMBER OF TIMES UNTIL THE ENERGY IS REDUCED BELOW 
+C	THE PUNCH-THROUGH ENERGY (PTE). IF THE INITIAL ENERGY WAS TOO SHORT
+C	THEN IT WON'T UNTIL PTE IS REACHED. IN THIS MOMENT, IPASS IS SET TO
+C	ONE, AND FROM THIS POINT EVERY FURTHER CALCULATION WILL IMPLY A
+C	REDUCTION BY HALF OF THE SIZE OF "DSEP", AND A CHANGE OF SIGN ONLY
+C	IF APROPRIATE.
+
+	IF (isign.LT.0) IPASS = 1		!IPASS=1 UNTIL "PTE" IS FOUND.
+	IF (IPASS.EQ.1) DSEP = isign * DSEP * 0.5
+
+	IF (ABS(DSEP).LT.0.05) THEN
+		EI = 0.00001
+		ISTAT = 0
+		GO TO 601
+	END IF
+
+	EP = EP + DSEP
+	ei = ep
+
+	GO TO 600
+
+1000	CONTINUE
+	GO TO 10
+2000	CONTINUE
+
+101	FORMAT(//////' PASSAGE OF CHARGED PARTICLE THROUGH ABSORBER',
+	1       ' SANDWICH '////'   AP = ',F6.0,'   ZP = ',F5.0,
+	1       '    INITIAL ENERGY = ',F12.5//)
+102	FORMAT('   ABSORBER SANDWICH CONTAINS - ',I2,' LAYERS'//)
+291	format(' start absorber calculations for   z =',f3.0
+     +,'  to z =',f3.0,/' and energies from emin =',f6.2
+     +,'  to emax =',f7.2,'  in ',f6.2,'MeV steps')
+311	FORMAT(//'  LAYER # ',I2,'    - SOLID ABSORBER -   ',
+	1' AREAL DENSITY = ',E10.4/' THICKNESS = ',E10.4,
+	1' CM   OR ',E10.4,' INCH      DENSITY =',E10.3,' G/CM3')
+312	FORMAT(//'   LAYER # ',I2,'    -  GAS  ABSORBER -   ',
+	1'  AREAL DENSITY = ',E10.4/'  PRESSURE = ',E10.4,
+	1' TORR    LENGTH ',E9.4,'CM    DENSITY =',E9.3,'MG/CM3')
+321	FORMAT(7X,' A =',F6.0,'  Z =',F5.0,'   AREAL DENSITY'
+	1,' (PARTIAL) = ',E12.5,' MG/CMSQ')
+401	FORMAT(' CALC IN-',I4,' STEPS'
+	1,'   ENERGY IN = ',F8.3,'    ENERGY OUT = ',F8.3
+	2,'(MEV)')
+402	FORMAT(' CHARGED PARTICLE STOPPED IN LAYER # ',I2)
+613	FORMAT (2X,'LAYER= ',I2,': ',I6,' ITERATIONS'
+	1,   ', E final= ',F10.4,' STATUS= ',I2)
+611	FORMAT (2X,'Ion  (A , Z): (',F4.0,' , ',F3.0
+	1,'), E(MeV)= ',F7.2,'  STOPPED IN LAYER ',I2/ 
+	1'   Esum = ',f7.2,'    Esum-E1 = ',f7.2,
+	1'   Esum-E1-E2 = ',f7.2)
+703	format('  eptable (',i2,', ',i3,', ',i1,' ) =',f8.2)
+
+	return
+	END
+
+	SUBROUTINE ADS(I1,SIGN,XN1,EPS,A,Z,E,DEE,ISTAT)
+
+C	SUBROUTINE FOR ENERGY LOSS CALCULATIONS
+C	CALL DEDX FRO STOPPING POWER CALCULATIONS
+
+	COMMON ISG(19),INN(19)
+	1      ,DEN(19),THK(19),PRS(19),XLN(19),ARDEN(19)
+	1      ,ZNUMB(19,4),ANUMB(19,4),ELNUM(19,4),CONCN(19,4)
+	1         ,PRDEN(19,4),PRTHK(19,4)
+
+C	N1= NUMBER OF SUBDIVISIONS FOR INTEGRATION
+C		OF ENERGY LOSS
+
+1000	CONTINUE
+	EH = E
+	N1 = IFIX(XN1+0.001)
+	DEDNEXT = 0.
+	DO 1010 K=1,N1
+		J1 = INN(I1)
+		ISGW = ISG(I1)
+		I = I1
+		DO 1001 J = 1,J1
+			AX = ANUMB(I,J)
+			ZX = ZNUMB(I,J)
+			FX = PRTHK(I,J)/XN1
+			DENST = PRDEN(I,J)
+			VH = VEL(EH,A)
+			CALL DEDX(Z,A,ZX,AX,DENST,EH,VH,ISGW,DEX,DE)
+			EH = EH + DE*SIGN*FX
+			IF(EH.LE.0.) THEN
+				IF (K.LE.2) THEN
+					N1 = N1 * 2
+					XN1 = FLOAT(N1)
+					GO TO 1000
+				END IF
+				ISTAT = -1
+				GO TO 9910
+			END IF
+			IF (K.LE.2) DEDNEXT = DEDNEXT + DE * FX
+1001		CONTINUE
+		IF (K.EQ.1) THEN
+			DED1ST = DEDNEXT
+			DEDNEXT = 0.0
+		END IF
+		IF (K.EQ.2) THEN
+			DDD = DED1ST - DEDNEXT
+			IF(DDD.LT.0.) DDD=-DDD
+			DDS = DED1ST + DEDNEXT
+			DDR = DDD/DDS
+			IF(DDR.GT.EPS) THEN
+				N1 = N1 * 2
+				XN1 = FLOAT(N1)
+				GO TO 1000
+			END IF
+		END IF
+1010	CONTINUE
+
+	ISTAT = 0
+9910	DEE = EH-E
+
+	RETURN
+	END
+
+	SUBROUTINE SETABS(INW,A,Z,AN,T,TH,D,DN)
+
+C	SUBROUTINE FOR SETTING UP COMPOSITE ABSORBER
+C	DATA (PARTIAL DENSITIES AND THICKNESSES)
+
+	DIMENSION A(4),Z(4),AN(4),T(4),D(4)
+
+	AW=0.
+	DO 1 I=1,INW
+	AW=AW+A(I)*AN(I)
+1	CONTINUE
+	DO 2 I=1,INW
+	AN(I)=A(I)*AN(I)/AW
+	T(I) = TH*AN(I)
+	D(I) = DN*AN(I)
+2	CONTINUE
+
+	RETURN
+	END
+
+	SUBROUTINE SETABG(INW,A,Z,AN,CN,T,TH,D,DN,PR,XL)
+
+C	SUBROUTINE FOR SETTING UP COMPOSITE ABSORBER DATA
+C	FOR GASEOUS LAYERS.
+
+	DIMENSION A(4),Z(4),AN(4),CN(4),T(4),D(4)
+	
+
+	P = PR/760.
+	X = XL/22.4
+
+	AWW=0.
+	AW=0.
+	DO 1 I=1,INW
+	AW = AW +A(I)*AN(I)
+	AWW= AWW+A(I)*AN(I)*CN(I)
+	T(I) = P*X*A(I)*AN(I)*CN(I)
+	D(I) = T(I)/XL
+1	CONTINUE
+	RETURN
+	END
+
+	FUNCTION VEL(ENER,A1)
+
+	VV=SQRT(2.13E-3*ENER/A1)
+	VEL=VV
+	RETURN
+	END
+
+	FUNCTION FKINEM(EP,AP,AT,TH)
+
+	IF(AP.GT.AT) GOTO 100
+	E=EP*AP**2/(AP+AT)**2
+	E=E*(COS(TH)+SQRT((AT/AP)**2-SIN(TH)**2))**2
+	FKINEM=E
+	RETURN
+100	FKINEM=0.
+	RETURN
+	END
+
+	FUNCTION FFKIN(EP,AP,AT,TH,Q)
+
+C	INELASTIC SCATTERING
+	B=AP**2*EP/(AP+AT)**2/(EP+Q)
+	D=AT**2/(AP+AT)**2*(1.+AP*Q/AT/(EP+Q))
+	E=(EP+Q)*B*(COS(TH)+SQRT(D/B-SIN(TH)**2))**2
+	FFKIN=E
+	RETURN
+	END
+
+	SUBROUTINE DEDX(Z1,A1,Z2,A2,RHO,ENER,V,IFG,DEDXHI,DEDXTO)
+
+C PROGRAM CALCULATES THE DIFFERENTIAL ENERGY LOSS DE/DX IN SOLID
+C TARGETS USING A SEMIEMPIRICAL FORMULA DEDUCED FROM EXPERIMENTAL
+C THE PROGRAM IS MODIFIED FOR GAS ABSORBERS.
+
+C REF.: K.BRAUNE,R.NOVOTNY,D.PELTE,D.HUSAR,D.SCHWALM,
+C PROCEEDINGS - SPRING MEETING OF THE GERMAN PHYSICAL
+C SOCIETY, VERHANDLUNGEN 4/1978
+
+C K.BRAUNE, DIPLOM, HEIDELBERG 1979
+
+
+C 	  H(Z2) IS A SUM OF FIVE GAUSSIAN FUNCTIONS.
+C 	       A1        MASS NUMBER  - PROJECTILE
+C 	       Z2        ATOMIC NUMBER ABSORBER
+C 	       A1        MASS NUMBER   ABSORBER
+C 	       RHO       DENSITY OF THE ABSORBER (GRAMM/CM**3)
+C 	                 (MEANLESS IF GAS ABSORBER )
+C 	       ENER      ENERGY OF THE PROJECTILE (MEV)
+C 	       V         VELOCITY OF THE PROJECTILE
+C 	                 IN MEV/(MG/CM**2)
+C       	Z1       ATOMIC NUMBER - PROJECTILE
+
+	IF(IFG.EQ.1) RHO=1.
+	XI=V**2/Z2
+C 	
+C	ABSORBER - FUNCTION
+C 	G(XI)=Y(EXP)-Y(THEORY) - IS DEDUCED FROM EXPERIMENTAL ENERGY LOSS
+C 	MEASUREMENTS.
+C 	
+C 	IF THE SAME ABSORBER WAS USED BEFORE , GO TO STATEMENT # 55
+
+	IF(A2.EQ.A2SAV.AND.Z2.EQ.Z2SAV) GOTO 55
+
+	Z2SAV=Z2
+	A2SAV=A2
+
+C  FUNCTION  Y
+
+	FY=54721.*(1.+5.15E-2*SQRT(A2/RHO)-EXP(-0.23*Z2))
+	IF(IFG.NE.1) GOTO 10
+	FY=54721.*(1.35-EXP(Z2*(-.13+.0014*Z2)))
+
+C  G(XI) IS THE DERIVATION OF A GASSIAN WITH VARIABLE HEIGHT H(Z2).
+
+10	IF(Z2.GT.26.) GOTO 20
+	G1=19.84*EXP(-.17*(Z2-4.25)**2)
+	GOTO 35
+20	G1=0.000001
+	IF(Z2.GT.38.) GOTO 40
+35	G2=17.12*EXP(-.12*(Z2-11.63)**2)
+	GOTO 50
+40	G2=0.0000001
+50	G3=7.95*EXP(-.015*(Z2-30.2)**2)
+	G4=5.84*EXP(-.022*(Z2-48.63)**2)
+	G5=7.27*EXP(-.005*(Z2-73.06)**2)
+	HZ2=(9.-(G1+G2+G3+G4+G5))*1.32E-5
+	ZWD=2./3.
+	Z2ZWD=Z2**ZWD
+
+C  MULTIPLICATIONFACTORS OF G(XI)
+
+	FG=1.2E-4*Z2*Z2+2.49E-2*A2/RHO
+	IF(IFG.NE.1) GOTO 52
+	FG=1.3/(1.+EXP(3.-Z2/5.))
+52	ALEFG=ALOG(2.7E-5/FG)
+
+C  CALCULATION OF G(XI)
+
+55	GXI=0.
+	IF(XI.GE.1.E-9.AND.XI.LE.5.E-4) THEN
+
+	SQXI=SQRT(XI)
+	C=2./Z2*SQXI/(1.+1.E4*SQXI)
+	IF(IFG.EQ.1) C=C/2.
+	FG0=1./(1.+(XI*10000.)**3)
+	AL=ALOG(XI)-ALEFG
+	GXI=(C-HZ2*AL*EXP(-.32*AL*AL))*FG0
+	ENDIF
+
+C  CALCULATION OF Y(XI)
+	Y=3.3E-4*ALOG(1.+XI*FY)+GXI
+C  ENERGY LOSS OF HEAVY IONS
+C  EFFECTIVE CHARGE
+	VV0=V*137.
+	FV=1.
+	IF(V.LE..62) FV=1.-EXP(-VV0)
+	AZ1=ALOG(1.035-.4*EXP(-.16*Z1))
+
+	QQ=V/Z1**.509
+	GHI=Z1
+	VZ1=(-116.79-3350.4*QQ)*QQ
+	IF(VZ1.GT.-85.2) GHI=Z1*(1.-EXP(VZ1))
+	IF(Z1.GT.2.) GHI=Z1*(1.-EXP(FV*AZ1-0.879*VV0/Z1**0.65))
+
+C  EFFECTIVE CHARGE FOR PROTONS AND ALPHA PARTICLES
+
+C ********************   RESULTS  ********************
+
+C  ELECTRONIC ENERGY LOSS DEDXHI
+
+	DEDXHI=GHI*GHI*Z2*Y/(A2*V**2)
+
+C  NUCLEAR ENERGY LOSS DEDXNU
+
+	ZA=SQRT(Z1**(ZWD)+Z2ZWD)
+	EPS=3.25E4*A2*ENER/(Z1*Z2*(A1+A2)*ZA)
+	SIGMAN=1.7*SQRT(EPS)*ALOG(EPS+2.1718282)/
+	1     (1.+6.8*EPS+3.4*EPS**1.5)
+	DEDXNU=SIGMAN*5.105*Z1*Z2*A1/(ZA*A2*(A1+A2))
+
+C  TOTAL ENERGY LOSS DEDXTO
+
+	DEDXTO=DEDXHI+DEDXNU
+
+	RETURN
+	END
+c

eloss_calculations/stopit/filtered.out

@@ -0,0 +1,6 @@
+Charged particle has Z =  2.  A =  4.  Initialenergy =    0.132 MeV.
+Energy lost in layer =   0.132 MeV
+Energy remaining =   0.000 MeV
+Charged particle has Z =  2.  A =  4.  Initialenergy =   10.132 MeV.
+Energy lost in layer =  10.133 MeV
+Energy remaining =   0.000 MeV

eloss_calculations/stopit/outfile.out

@@ -0,0 +1,50 @@
+4
4
+
+
+Charged particle has Z =  2.  A =  4.  Initialenergy =    0.132 MeV.
+
+
+Absorber contains  1 layers.
+
+
+Layer #  1 -  Solid Absorber
+Density = 0.2329E+01 g/cm3      Thickness = 0.2329E+02      mg/cm2
+Z = 14.  A = 28.    Atoms per molecule =  1.
+Energy lost in layer =   0.132 MeV
+
+Energy remaining =   0.000 MeV
+
+Output also written to stopit.log
+
+  1 define stopee
+  2 define absorber
+  3 edit absorber
+  4 run with current parameters
+  5 find thickness of absorber to stop the stopee
+  6 print status of data
+  7 stop
+4
4
+
+
+Charged particle has Z =  2.  A =  4.  Initialenergy =   10.132 MeV.
+
+
+Absorber contains  1 layers.
+
+
+Layer #  1 -  Solid Absorber
+Density = 0.2329E+01 g/cm3      Thickness = 0.2329E+02      mg/cm2
+Z = 14.  A = 28.    Atoms per molecule =  1.
+Energy lost in layer =  10.133 MeV
+
+Energy remaining =   0.000 MeV
+
+Output also written to stopit.log
+
+  1 define stopee
+  2 define absorber
+  3 edit absorber
+  4 run with current parameters
+  5 find thickness of absorber to stop the stopee
+  6 print status of data
+  7 stop

eloss_calculations/stopit/script.exp

@@ -0,0 +1,106 @@
+#!/usr/bin/expect -f
+#
+set timeout -1
+spawn ./a.out
+match_max 100000
+expect -exact "\r
+  1 define stopee\r
+  2 define absorber\r
+  3 edit absorber\r
+  4 run with current parameters\r
+  5 find thickness of absorber to stop the stopee\r
+  6 print status of data\r
+  7 stop\r
+"
+send -- "1\r"
+expect -exact "1\r
+ Enter Z and A of stopee.\r
+"
+send -- "2\r"
+expect -exact "2\r
+"
+send -- "4\r"
+expect -exact "4\r
+ Enter energy in MeV.\r
+"
+send -- "5.486\r"
+expect -exact "5.486\r
+\r
+  1 define stopee\r
+  2 define absorber\r
+  3 edit absorber\r
+  4 run with current parameters\r
+  5 find thickness of absorber to stop the stopee\r
+  6 print status of data\r
+  7 stop\r
+"
+send -- "2\r"
+expect -exact "2\r
+\r
+ How many stopping layers are there in the absorber?\r
+"
+send -- "1\r"
+expect -exact "1\r
+\r
+Is layer  1 a solid, gas, or standardized medium?\r
+(0 = solid, 1 = gas, 2 = standardized medium)\r
+(Standardized media are: CO2, Si, C \[graphite\], C4H10, CF4,\r
+CH2, CD2, He, and H2)\r
+"
+send -- "2\r"
+expect -exact "2\r
+Which standardized medium for layer  1?\r
+1: CO2\r
+2: Si\r
+3: C (graphite)\r
+4: C4H10\r
+5: CF4\r
+6: CH2\r
+7: CD2\r
+8: He-gas\r
+9: H2-gas\r
+"
+send -- "2\r"
+expect -exact "2\r
+Enter thickness(microns) for layer  1\r
+"
+send -- "100\r"
+expect -exact "100\r
+\r
+  1 define stopee\r
+  2 define absorber\r
+  3 edit absorber\r
+  4 run with current parameters\r
+  5 find thickness of absorber to stop the stopee\r
+  6 print status of data\r
+  7 stop\r
+"
+send -- "4\r"
+expect -exact "4\r
+\r
+\r
+Charged particle has Z =  2.  A =  4.  Initialenergy =    5.486 MeV.\r
+\r
+\r
+Absorber contains  1 layers.\r
+\r
+\r
+Layer #  1 -  Solid Absorber\r
+Density = 0.2329E+01 g/cm3      Thickness = 0.2329E+02      mg/cm2\r
+Z = 14.  A = 28.    Atoms per molecule =  1.\r
+Energy lost in layer =   5.486 MeV\r
+\r
+Energy remaining =   0.000 MeV\r
+\r
+Output also written to stopit.log\r
+\r
+  1 define stopee\r
+  2 define absorber\r
+  3 edit absorber\r
+  4 run with current parameters\r
+  5 find thickness of absorber to stop the stopee\r
+  6 print status of data\r
+  7 stop\r
+"
+send -- "7\r"
+expect eof

eloss_calculations/stopit/stopit.expect

@@ -0,0 +1,106 @@
+#!/usr/bin/expect -f
+#
+# This Expect script was generated by autoexpect on Sat Jan 31 17:06:49 2026
+# Expect and autoexpect were both written by Don Libes, NIST.
+#
+# Note that autoexpect does not guarantee a working script.  It
+# necessarily has to guess about certain things.  Two reasons a script
+# might fail are:
+#
+# 1) timing - A surprising number of programs (rn, ksh, zsh, telnet,
+# etc.) and devices discard or ignore keystrokes that arrive "too
+# quickly" after prompts.  If you find your new script hanging up at
+# one spot, try adding a short sleep just before the previous send.
+# Setting "force_conservative" to 1 (see below) makes Expect do this
+# automatically - pausing briefly before sending each character.  This
+# pacifies every program I know of.  The -c flag makes the script do
+# this in the first place.  The -C flag allows you to define a
+# character to toggle this mode off and on.
+
+set force_conservative 0  ;# set to 1 to force conservative mode even if
+			  ;# script wasn't run conservatively originally
+if {$force_conservative} {
+	set send_slow {1 .1}
+	proc send {ignore arg} {
+		sleep .1
+		exp_send -s -- $arg
+	}
+}
+
+#
+# 2) differing output - Some programs produce different output each time
+# they run.  The "date" command is an obvious example.  Another is
+# ftp, if it produces throughput statistics at the end of a file
+# transfer.  If this causes a problem, delete these patterns or replace
+# them with wildcards.  An alternative is to use the -p flag (for
+# "prompt") which makes Expect only look for the last line of output
+# (i.e., the prompt).  The -P flag allows you to define a character to
+# toggle this mode off and on.
+#
+# Read the man page for more info.
+#
+# -Don
+
+
+set timeout -1
+spawn ./a.out
+match_max 100000
+expect "*  7 stop\r
+"
+send -- "1\r"
+expect -exact "1\r
+ Enter Z and A of stopee.\r
+"
+send -- "4"
+expect -exact " "
+send -- "2\r"
+expect -exact "2\r
+"
+send -- "4\r"
+expect -exact "4\r
+ Enter energy in MeV.\r
+"
+send -- "5.486\r"
+expect -exact "5.486\r
+\r
+  1 define stopee\r
+  2 define absorber\r
+  3 edit absorber\r
+  4 run with current parameters\r
+  5 find thickness of absorber to stop the stopee\r
+  6 print status of data\r
+  7 stop\r
+"
+send -- "2\r"
+
+expect "*How many stopping layers are there in the absorber?\r"
+send -- "1\r"
+
+expect  "*CH2, CD2, He, and H2)\r"
+send -- "2\r"
+
+expect  "*Which standardized medium for layer  1?\r
+1: CO2\r
+2: Si\r
+3: C (graphite)\r
+4: C4H10\r
+5: CF4\r
+6: CH2\r
+7: CD2\r
+8: He-gas\r
+9: H2-gas\r
+"
+send -- "2\r"
+expect  "
+* Enter thickness(microns) for layer  1\r
+"
+send -- "100\r"
+expect  "
+*  7 stop\r
+"
+send -- "4\r"
+expect  "
+*  7 stop\r
+"
+send -- "7\r"
+expect eof

eloss_calculations/stopit/stopit_pex.py

@@ -0,0 +1,87 @@
+import pexpect as px
+import sys,time,os
+import matplotlib.pyplot as plt
+
+if len(sys.argv) != 8:
+    print("Usage: python3 stopit_pex.py <z> <a> <thickness_um_layer1> <thickness_um_layer2> <energy_start_MeV>  <energy_stop_MeV> <energy_step>")
+    print("Modify the \#absorber section in script to change stopping material")
+    quit()
+
+stoppeeZ=int(sys.argv[1])
+stoppeeA=int(sys.argv[2])
+layer1_thickness = float(sys.argv[3]) #um
+layer2_thickness = float(sys.argv[4]) #um
+min_stoppee_energy=float(sys.argv[5])
+max_stoppee_energy=float(sys.argv[6])
+step_stoppee_energy=float(sys.argv[7])
+
+#plotting prep
+'''
+plt.ion()
+fig=plt.figure()
+ax=fig.add_subplot(111)
+ax.autoscale_view(True,True,True)
+line,=ax.plot([],[],'b-')
+ax.set_xlim(0,0.5*max_stoppee_energy)
+ax.set_ylim(0,0.05*max_stoppee_energy)
+'''
+stopit_child = px.spawn("./a.out")
+stopit_child.logfile = None
+
+def run_stopit(child,Z,A,Energy_MeV,thickness_um):
+	child.expect("7 stop\r")
+
+	#stoppee
+	child.send("1\r")
+	child.expect("Enter Z and A of stopee.\r")
+	child.send(str(Z)+" "+str(A)+"\r")
+	child.expect("Enter energy in MeV.\r")
+	child.send(str(Energy_MeV)+"\r")
+
+	#absorber - standard medium, Si
+	child.expect("7 stop\r")
+	child.send("2\r")
+	child.expect("\r")
+	child.send("1\r")
+	child.expect("He, and H2\)\r")
+	child.send("2\r")
+	child.expect("H2-gas\r")
+	child.send("2\r") #first 2 for 'define absorber', then 1 for '1 layer', then 2 for 'standardized material', then '2' for Si.
+	child.expect("\r")
+	child.send(str(thickness_um)+"\r")
+
+	#run!
+	child.expect("7 stop\r")
+	child.send("4\r")
+	child.expect("stopit.log\r")
+	#get all the output that shows up before the string 'stopit.log' above.
+	#the split() splits it up at all whitespaces. we can pick out the substrings of interest from this bunch now
+	results = child.before.decode('utf-8').split()
+
+	return [results[12],results[48],results[53]]
+
+layer1_edep = []
+layer2_edep = []
+stoppee_energy = min_stoppee_energy #in MeV
+while(stoppee_energy < max_stoppee_energy):
+	[in_en,e_lost ,e_left ] = run_stopit(stopit_child,stoppeeZ,stoppeeA,stoppee_energy,layer1_thickness)
+	[in_en,e_lost2,e_left2] = run_stopit(stopit_child,stoppeeZ,stoppeeA,e_left        ,layer2_thickness)
+
+	print(stoppee_energy,e_lost,e_lost2)
+	sys.stdout.flush()
+	layer1_edep.append(float(e_lost))
+	layer2_edep.append(float(e_lost2))
+	'''
+	line.set_data(layer2_edep,layer1_edep) #dE in y-axis, E in x-axis when possible
+	ax.relim()
+	plt.draw()
+	plt.pause(0.004)
+	'''
+	stoppee_energy += step_stoppee_energy
+
+stopit_child.expect("7 stop\r")
+stopit_child.send("7\r")
+stopit_child.expect(px.EOF)
+stopit_child.close()
+
+plt.show(block=True)

eloss_calculations/test_tgraph.C

@@ -0,0 +1,17 @@
+void test_tgraph() {
+    TGraph MeV_to_cm("alphas_in_250torr_mix_filtered.txt","%lf %*lf %lf");
+    TGraph cm_to_MeV(MeV_to_cm.GetN(), MeV_to_cm.GetY(), MeV_to_cm.GetX());
+
+    MeV_to_cm.Draw("AL");
+    while(gPad->WaitPrimitive());
+    cm_to_MeV.Draw("AL");
+    gPad->Modified(); gPad->Update();
+    while(gPad->WaitPrimitive());
+
+    double e=6.0,L=4.0;
+    while(e>5.0) {
+        double temp = MeV_to_cm.Eval(e)-L;
+        std::cout << "If we detect an alpha at " << e << " MeV, covering a path of " << L << " cm, " << cm_to_MeV.Eval(MeV_to_cm.Eval(e)-L) << " MeV is its original energy " << L << " cm prior" << std::endl;
+        e-=0.05;
+    }
+} 

gmsx3/func1.h

@@ -0,0 +1,41 @@
+#include "TF1.h"
+
+double model2(double *x, double *par) {
+    /*  'Potential Well' of width 2a from from xx-a to xx+a
+        xx is coordinate about the point of origin, set at x=center
+        v0 is the y-offset of the potential
+        k is the 'steepness' of the potential
+
+        continuous across xx-a and xx+a, and differentiable
+    */
+
+    double center= par[3];
+    double xx = x[0]-center;
+    double a = TMath::Abs(par[0]);
+    double k = TMath::Abs(par[1]);
+    double v0 = par[2];
+
+    if(xx < -a)
+        return k*(xx+a)*(xx+a) + v0;
+    else if(xx > a)
+        return k*(xx-a)*(xx-a) + v0;
+    else
+        return v0;
+}
+
+void func1() {
+    //TF1 f1("bowl",model,-2.,2.,2);
+    TCanvas c("c1","c1",800,600);
+    TF1 f1("bowl",model2,-10.,10.,4);
+    f1.SetMaximum(10);
+
+    for(int i=-4; i<4; i++) {
+        f1.SetParameters(.4,100,2,i); //a, k, v0, center
+        f1.SetNpx(100000);
+        if(i==-4) f1.Draw("L");
+        f1.DrawCopy("L SAME");
+        c.Modified(); c.Update();
+        //c.SaveAs(Form("%d.png",out));
+        while(c.WaitPrimitive());
+   }
+}

gmsx3/intgm_sx3.h

@@ -0,0 +1,216 @@
+#include "../Armory/HistPlotter.h"
+#include <Minuit2/FCNBase.h>
+#include <Math/Minimizer.h>
+#include <Math/Factory.h>
+#include <Math/Functor.h>
+#include <TMath.h>
+#include <TPad.h>
+#include <cassert>
+#include <vector>
+#include <array>
+#include <iostream>
+#include <TF1.h>
+#include "func1.h"
+static long iters=0;
+
+//class intgm_sx3 : public ROOT::Minuit2::FCNBase {
+class intgm_sx3 {
+        int N;
+
+        //L.at(0).at(3).at(n) is front strip = 0, back pad = 3, nth datapoint
+        std::array<std::array<std::vector<double>,4>,4> L,R,B;
+        //std::array<std::array<double,5>,4> stripedge; //stripedge.at(i).at(j) is the jth edge of the ith strip. there are five edges for the four strips 'i'=0 to 3, (0,1) (1,2) (2,3) (3,4) for each
+        //the edges are at  -2a, -a, 0, a, 2a respectively so we enforce four ratios in the chi2 value - 'a' can be held constant, no need to fit it.
+
+        //assume z = M*(aL-bR)
+        //stripedge[i][1] = max(z) when pad==0 = min(z) when pad==1  this should be -1
+        //stripedge[i][2] = max(z) when pad==1 = min(z) when pad==2. this should be 0
+        //stripedge[i][3] = max(z) when pad==2 = min(z) when pad==3. this should be 1
+
+        //i.e. stripedge[i][j] = max(z) when pad == j-1, min(z) when pad==j, for i= 1,2,3
+
+        //ncounts.at(frontch).at(backch) is the number of (L,R,B) tuples we've filled (frontch,backch) coordinates in the detector
+        std::array<std::array<long,4>,4> ncounts;
+        TH1F *localhists[4][4]; //one histogram for each fc, bc combination
+        HistPlotter *plotter;
+        TF1 *pos_weight[4];
+        TF1 *energywell;
+    public:
+        intgm_sx3() {
+            for(int bc=0; bc<4; bc++) {
+                for(int fc=0; fc<4; fc++) {
+                    L[fc][bc].reserve(1000);
+                    R[fc][bc].reserve(1000);
+                    B[fc][bc].reserve(1000);
+                    //localhists[fc][bc] = new TH1F(Form("h_%d_%d",fc,bc),Form("h_%d_%d",fc,bc),1000,-4.,4.);
+                    ncounts[fc][bc] = 0;
+                }
+                pos_weight[bc] = new TF1(Form("b_strip%d",bc),model2,-10,10,4); //from -10, to 10, 4 parameters
+                pos_weight[bc]->SetParameters(1.0,10,1.,3-2*bc); //centers at 1, 3.,5,7 Width 2a with a=1.0
+                pos_weight[bc]->SetNpx(1'000'000);
+            }
+            energywell = new TF1("ewell",model2,0,2000,4); //0 to 2000 channels, 4 params
+            energywell->SetParameters(1000,20,1,1500); //center the back E values at 1500 +/- 500
+            energywell->SetNpx(1'000'000);
+            N=0;
+        }
+        void set_plotter(HistPlotter *p) {plotter=p;}
+        void set_iters(long i) { iters=i;}
+        intgm_sx3(HistPlotter *p) : plotter(p) {
+            for(int bc=0; bc<4; bc++) {
+                for(int fc=0; fc<4; fc++) {
+                    L[fc][bc].reserve(1000);
+                    R[fc][bc].reserve(1000);
+                    B[fc][bc].reserve(1000);
+                    //localhists[fc][bc] = new TH1F(Form("h_%d_%d",fc,bc),Form("h_%d_%d",fc,bc),1000,-4.,4.);
+                    ncounts[fc][bc] = 0;
+                }
+                pos_weight[bc] = new TF1(Form("b_strip%d",bc),model2,-10,10,4); //from -10, to 10, 4 parameters
+                                                //a/2, k, v0, center
+                pos_weight[bc]->SetParameters(0.92,10,1.,-1.*(3-2*bc)); //centers at 7, 5.,3,1 Width 2a with a=1.0
+                pos_weight[bc]->SetNpx(1'000'000);
+            }
+            energywell = new TF1("ewell",model2,0,8000,4); //0 to 2000 channels, 4 params
+//            energywell->SetParameters(60,10,0,1430); //center the back E values at 1430 +/- 60
+            energywell->SetParameters(400,10,0,5246); //center the back E values at 5486 +/- 600
+            energywell->SetNpx(1'000'000);
+            N=0;
+        }
+
+        inline void fill(int fc, int bc, double leftE, double rightE, double backE) {
+            /*
+             *
+             */
+            assert(fc>=0 && fc<=3 && "Front channels should fit the range 0 to 3 inclusive!");
+            assert(bc>=0 && bc<=3 && "Back channels should fit the range 0 to 3 inclusive!");
+            if(leftE>0 && rightE >0 && backE>0) {
+                L[fc][bc].emplace_back(leftE);
+                R[fc][bc].emplace_back(rightE);
+                B[fc][bc].emplace_back(backE);
+                ncounts[fc][bc]+=1;
+                N+=1;
+            }
+        }
+
+        inline void print() {
+            for(int i=0; i<16; i++) {
+                std::cout << ncounts[i%4][i/4] << std::endl;
+            }
+        }
+        inline void plot(std::string comment, const double* params) {
+            std::array<double,4> l,r,b,bo,ro,lo,offset,stretch; //aliases to help with book-keeping
+            std::array<std::array<double,4>,4> back_gains;// back_gains[fc][bc] are for fc,bc firing in combo
+            for(int ctr=0; ctr<4; ctr++) {
+                r[ctr] = params[ctr];
+            }
+            for(int ctr=4; ctr<20; ctr++) {
+                int bch = (ctr-4)%4;
+                int fch = (ctr-4)/4;
+                back_gains[bch][fch] = params[ctr];
+            }
+            for(int ctr=20; ctr<24; ctr++) {
+                stretch[ctr-20] = params[ctr];
+            }
+            for(int ctr=24; ctr<28; ctr++) {
+                l[ctr-24] = params[ctr];
+            }
+            for(int fc=0; fc<4; fc++) {
+                for(int bc=0; bc<4; bc++) {
+                    for(int n=0; n<ncounts[fc][bc]; n++) {
+                        if(plotter) {
+                            double left = l[fc]*L[fc][bc].at(n);
+                            double right = r[fc]*R[fc][bc].at(n);
+                            double back = back_gains[bc][fc]*B[fc][bc].at(n);
+                            //double zpos = (left - right)/(left+right);
+                            double zpos = stretch[fc]*(left - right)/(left+right);// + offset[fc]; //back;
+                            plotter->Fill2D(Form("normlf_fc%d_%d_%s",fc,bc,comment.c_str()),800,0,1.,800,0,1.,left/back, right/back,"l_vs_r");
+                            plotter->Fill2D(Form("normlf_all_%s",comment.c_str()),800, 0, 1., 800, 0, 1.,left/back, right/back);
+                            plotter->Fill2D(Form("case_f%d_b%d_%s",fc,bc,comment.c_str()),800,0,8192,800,0,8192,left+right,back,"l_vs_r");
+                            plotter->Fill2D(Form("case_all_%s",comment.c_str()),800,0,8192,800,0,8192,left+right,back);
+                            //plotter->Fill2D(Form("z_vs_backe_f%d_b%d_%s",fc,bc,comment.c_str()),800,-10,10,800,0,8192,zpos,back,"z_vs_be");
+                            plotter->Fill2D(Form("z_vs_backe_all_%s",comment.c_str()),800,-10,10,800,0,8192,zpos,back);
+                        } //end if plotter
+                    }// end for-n
+                }//end for-bc
+            }//end for-fc
+        }//end plot()
+
+    //    double operator()(const std::vector<double>& params) const override{
+        double eval(const double* params) const {
+            iters+=1;
+
+            std::array<double,4> l,r,b,bo,ro,lo, offset, stretch; //aliases to help with book-keeping
+            std::array<std::array<double,4>,4> back_gains;// back_gains[fc][bc] are for fc,bc firing in combo
+            for(int ctr=0; ctr<16; ctr++) {
+                int bch = (ctr)%4;
+                int fch = (ctr)/4;
+                back_gains[bch][fch] = params[ctr];
+            }
+            for(int ctr=16; ctr<20; ctr++) {
+                r[ctr-16] = params[ctr];
+                l[ctr-16] = 1.0;
+            }
+            for(int ctr=20; ctr<24; ctr++) {
+            	stretch[ctr-20] = params[ctr];
+            }
+            double result=0, sumcount=0;
+            for(int fc=0; fc<4; fc++) {
+                for(int bc=0; bc<4; bc++) {
+                    //if(bc >= 1 || fc >= 1 ) continue;
+                    if(ncounts[fc][bc] == 0 && iters ==0) {
+                        std::cout << "Missing any data in front:" << fc << " back:" << bc << " combination." << std::endl;
+                    }
+
+                    for(int n=0; n<ncounts[fc][bc] ; n++) {
+                        //double left = l[fc]*L[fc][bc].at(n) + lo[fc];
+                        //double right = r[fc]*R[fc][bc].at(n) + ro[fc];
+                        //double back = b[bc]*B[fc][bc].at(n) + bo[bc];
+                        //double add = TMath::Power(left + right - back,2);
+                        if(B[fc][bc].at(n)<100) continue;//ignore events too close to noise threshold
+
+                        double left = l[fc]*L[fc][bc].at(n);
+                        double right = r[fc]*R[fc][bc].at(n);
+                        double back = back_gains[bc][fc]*B[fc][bc].at(n);
+                        double lnorm = left/B[fc][bc].at(n);
+                        double rnorm = right/B[fc][bc].at(n);
+
+                        //double add = TMath::Power(left/back + right/back - 1.0,2);
+                        double add = TMath::Power(left + right - back,2);
+                        double zpos = stretch[fc]*(left - right)/(left+right); //back;
+                        std::cout << zpos << " " << pos_weight[bc]->Eval(zpos) << " " << bc << std::endl;
+                        double add_position = pos_weight[bc]->Eval(zpos);
+                        double eback_align_penalty = energywell->Eval(back);
+/*                      if(back>1000) zmid[fc][bc] += zpos;
+                        if(back>1000 && zpos < zmin[fc][bc]) zmin[fc][bc] = zpos;
+                        if(back> 1000 && zpos > zmax[fc][bc]) zmax[fc][bc] = zpos;
+
+                        if(back>1000) {
+                            localhists[fc][bc]->Fill(zpos);
+                        }*/
+
+                        result += add_position;
+                        //result += add;
+                        result += eback_align_penalty;
+                        sumcount+=1;
+                        //if(bc==0) std::cout << add << " " << add_position << " " <<  zpos << std::endl;
+                        //To avoid drift towards (0,0,0) trivial solution. This value ~1 close to (1,1,1)
+                        //result+=(1e-3/(TMath::Power(l[fc],2)+TMath::Power(r[fc],2)+TMath::Power(b[bc],2)+1e-9));
+						//result+=(1e-3/(TMath::Power(l[fc],2)+TMath::Power(r[fc],2)+TMath::Power(b[bc],2)+1e-9));
+                    } //end for-n
+                } //end for-bc
+            } //end for-fc
+            result/=sumcount; //normalize, so the value doesn't scream
+            if(iters%1'000==0) {
+                std::cout << "iters : " << iters << " params: " << std::endl ;
+                for(int i=0 ; i< 10; i++) std::cout << params[i] << " " << std::flush;
+                std::cout<< std::endl;
+                for(int i=10 ; i< 20; i++) std::cout << params[i] << " " << std::flush;
+                std::cout << std::endl <<  " result: " << result << std::endl;
+            } //end if
+            return result;
+        } //end eval()
+
+        //double Up() const override { return 1.0; } // Required by minuit2 FCBase
+
+};
+

gmsx3/utilities_orr.h

@@ -0,0 +1,338 @@
+#ifndef UTILS_ORR_H
+#define UTILS_ORR_H
+#include "datatypes.h"
+#include "HistPlotter.h"
+#include "Geometry_orr.h" //contains orruba geometry constants
+#include <cassert>
+#include <stdio.h>
+#include <cassert>
+#include <cstdint>
+#include <fcntl.h>
+#include <unistd.h>
+#include <vector>
+//#include <cstring>
+#include <string>
+#include <fstream>
+#include <sstream>
+#include <iostream>
+#include <cmath>
+#include <iomanip>
+#include <set>
+#include <sys/stat.h>
+#include "TMath.h"
+#include "counters.h"
+
+static named_counter oc("named orruba counters");
+
+inline float get_filesize(std::string filename) {
+    struct stat st;
+    stat(filename.c_str(), &st);
+    return st.st_size;
+}
+
+class orruba_params {
+    /*
+    *   
+    */
+public:
+    int chnum=DEFAULT_NULL; //!< global channel number
+    std::string type = "-";
+    int id=DEFAULT_NULL;
+    int layer=DEFAULT_NULL;
+    int frontback=DEFAULT_NULL;
+    int updown=DEFAULT_NULL;
+    int subid=DEFAULT_NULL;
+    int leftright=DEFAULT_NULL;
+
+    float ped=DEFAULT_NULL;
+    float offset=DEFAULT_NULL;
+    float gain=DEFAULT_NULL;
+    float gain2=DEFAULT_NULL; //for use with sx3's
+};
+
+class sx3_geometry_scalefactors {
+public:
+    //If sx3 has L, R being the left and right extremities, we choose add, stretch here such that
+    // x_mm = (x_raw+add)*stretch; so add=abs(L), stretch=75/(abs(L)+R)
+    float add[4];
+    float stretch[4];
+};
+
+class qqq5_finegains {
+public:
+    std::array<std::pair<float,float>,32> front;
+    //front.at(30).first = slope at clkpos 0, ring 30 for E front layer
+    //front.at(30).second = intercept for the same as above
+    std::array<std::pair<float,float>,4> back;
+};
+class sx3_fbgains {
+public:
+    //Order of indices is [pad][strip]
+    float padoffsets[4][4];
+    float padgains[4][4];
+
+    float stripLoffsets[4][4];
+    float stripLgains[4][4];
+
+    float stripRoffsets[4][4];
+    float stripRgains[4][4];
+};
+
+//Metadata ORRUBA needs to know about itself, to be configured at the start
+extern std::array<orruba_params,MAX_ORRUBA_CHANS> o_params;
+extern std::array<sx3_fbgains,24> sx3_xtalk_gains; //every sx3 needs to be gainmatched as a frontL-back, frontR-back pair (pad strip pair)
+extern std::array<sx3_geometry_scalefactors,24> sx3gs;
+extern std::array<qqq5_finegains,4> qqq5_fg_dE, qqq5_fg_E;
+
+class type19Raw {
+public:
+    long int timestamp;
+    std::vector<unsigned short> ch;
+    std::vector<unsigned short> val;
+
+    type19Raw() : timestamp(0), ch(50,0), val(50,0) {} //Reserve 50 for size of these vectors, initial value of zero
+    void print() const {
+        /*
+            print()
+            Prints type19Raw's contents to stdout for monitoring
+        */
+        std::cout << "------" << std::endl;
+        for(unsigned int ii=0; ii<ch.size(); ii++) {
+            std::cout << ch.at(ii) << " " << val.at(ii) << std::endl;
+        }
+    }
+};
+
+class sx3 {
+public:
+    //TODO: Convert to std::array
+    //Holds all information in an event, including ped subtraction+scaling. back[2].at(0) will have the largest energy seen in ch2, if any
+    std::vector<float> back[4];
+    std::vector<float> frontL[4];
+    std::vector<float> frontR[4];
+
+    double ts = DEFAULT_NULL;
+    //Easy lookup of final calibrated event. Only filled for valid cases, assumed for now to be 1L, 1R, 1B
+    float frontX=DEFAULT_NULL;
+    float frontXmm=DEFAULT_NULL;
+    float frontE=DEFAULT_NULL;
+    float backE=DEFAULT_NULL;
+    int stripF=DEFAULT_NULL;
+    int stripB=DEFAULT_NULL;
+    float frontEL=DEFAULT_NULL;
+    float frontER=DEFAULT_NULL;
+
+    float phi=DEFAULT_NULL; //
+
+    std::set<int> valid_front_chans;
+    std::set<int> valid_back_chans;
+    std::set<int> unmatched_front_chans; //every front channel is unmatched and invalid at first. when it gets matched, it gets removed and sent to valid
+
+    bool foundevent=false;
+    bool valid=false;//valid will be set to false in all cases where we have ambiguity
+    int flags=-1;//flags settable to different types of values to indicate different invalid situations
+
+    void fillevent(const std::string& position, const int subchannel, const float value); //make 'const' what functions don't need to change, helps with performance
+    void validate(const sx3_fbgains&, const sx3_geometry_scalefactors&);
+};
+
+class qqq5 {
+public:
+    //Holds all information in an event, including ped subtraction+scaling. front[2].at(0) will have the largest energy seen in ch2, if any
+    //TODO: Convert to std::array
+    std::vector<float> back[4];
+    std::vector<float> front[32];
+
+    double ts = DEFAULT_NULL;
+    float selftheta=DEFAULT_NULL,selfrho=DEFAULT_NULL;
+    //Easy lookup of the final calibrated event. Only filled for valid cases.
+    double frontE=DEFAULT_NULL;
+    double backE=DEFAULT_NULL;
+    int frontch;
+    int backch;
+
+    std::pair<int,int> adj_front_strips = {-1,-1};
+    std::pair<int,int> adj_back_strips = {-1,-1};
+    std::set<int> valid_front_chans; //list of channels that fire. can inspect size() of these to see if there are front/back events
+    std::set<int> valid_back_chans; // we use std::set since it makes for very readable code
+
+    bool foundevent=false;
+    bool valid=false; //valid will be set to false in all cases where we have ambiguity
+    int flags=-1; //flags settable to different types of values to indicate different invalid situations
+
+    void fillevent(const std::string& position, const int subchannel, const float value); //make 'const' what functions don't need to change, helps with performance
+    void validate();
+};
+
+struct orrubaevent {
+    //Every clean, valid charged-particle event will have these four parts
+    float dE=DEFAULT_NULL;  //!< true energy-loss in the dE layer. Found by gainmatching ADC readout to alpha data
+    float E=DEFAULT_NULL;   //!< energy deposited in the  E layer. When summed with dE, gives true energy in keV deposited by the particle in ORRUBA
+    float dE_PID = DEFAULT_NULL; //!< dE scaled for dE-layer's thickness, reducing the spread due to angular straggling by explicitly accounting for it. This will give a sharper pid plot that can be gated on better
+    float dE_linPID = DEFAULT_NULL; //!< dE_PID, but linearized using the prescription described in, say, PhysRevC.90.034601 (2014). dE_linPID = ((dE+E)^a-E^a)^(1/a) where a ~ 1.68 is chosen empirically
+    float Theta=DEFAULT_NULL; //!< Laboratory polar angle of event in radians, deprecated
+    float Phi=DEFAULT_NULL; //!< Lab azimuthal angle of event in radians, deprecated
+
+    //Helpful indices to make dE-E plots
+    std::string type; //!< "endcap" vs "barrel"
+    //!< Identify the position of the detector in the barrel, usually in accordance with the channel map: say we might learn detector is at "Quad 4" or "clk_pos 10", together with 'type'. Useful with HistPlotter class
+    int position=DEFAULT_NULL;
+	int subchdE_1=DEFAULT_NULL, subchdE_2=DEFAULT_NULL; //!< Identify the subchannels corresponding to the two sides of the dE detector. To avoid confusion, 1=strip(sx3), ring(qqq) and 2=pad(sx3), wedge(qqq)
+	int subchE_1=DEFAULT_NULL, subchE_2=DEFAULT_NULL; //!< Identify the subchannels corresponding to the two sides of the E detector. Same convention as above
+
+    float x=DEFAULT_NULL,y=DEFAULT_NULL,z=DEFAULT_NULL; //!< Laboratory x,y,z coordinates of the event from the E layer
+    float r0=DEFAULT_NULL,theta0=DEFAULT_NULL,phi0=DEFAULT_NULL; //!< vector elements from hit to origin for E layer
+    float r1=DEFAULT_NULL,theta1=DEFAULT_NULL,phi1=DEFAULT_NULL; //!< vector elements from hit to origin for dE layer
+};
+
+/*TODO:
+    * There will be some use for a class such that it stores 
+        PhysicsEvent = <Ex, Brho, THeta4, ...orrubaevent >
+    * Once the 'orrubaevent' structs are made, it should be callable. like 
+    std::vector<PhysicsEvent> getPhysicsFromVertices(Kinematics dpkin, std::vector<orrubaevent> orvec);?
+    * should the physics just go sit in orrubaevent? maybe orruba2024 class can have a kinematics object in it?
+*/
+
+class orruba2024 {
+private:
+//Class expected to be changed for each version of the analysis code
+public:
+    bool found_trk, found_trkpresc, found_tdcq, found_s800e1, found_s800trg,
+         found_rf, found_gt, found_si, found_siup;
+
+	bool found_de, found_e, found_qqq, found_sx3;//orruba
+    long long timestamp=DEFAULT_NULL;
+    std::vector<type19Raw> o_rawvec;
+
+    std::array<qqq5,4> uendcapE;
+    std::array<qqq5,4> uendcapdE;
+    std::array<sx3,12> ubarrelE;
+    std::array<sx3,12> ubarreldE;
+
+    std::array<sx3,2> dbarrelE;
+    //Results after post-processing, including possible multiplicities
+    std::vector<orrubaevent> events;
+
+    double target_z_offset;
+
+    //void initialize(const std::string & filename1, const std::string& filename2, const std::string& filename3);
+	float tdc_trk = DEFAULT_NULL;
+	float tdc_trksi = DEFAULT_NULL; //trackerpr+si stops are combined into one channel
+	float tdc_trkp = DEFAULT_NULL;
+	float tdc_q = DEFAULT_NULL;
+	float tdc_s800e1 = DEFAULT_NULL;
+	float tdc_s800trg = DEFAULT_NULL;
+	float tdc_rf = DEFAULT_NULL;
+	float tdc_rf_unwrap = DEFAULT_NULL;
+	float tdc_gt = DEFAULT_NULL;
+	float tdc_si = DEFAULT_NULL;
+	float tdc_siup = DEFAULT_NULL;
+    orruba2024(const std::vector<type19Raw>& orvec);
+    void postprocess(); //generates 'events', and performs validations on freshly entered data. performs fine gainmatching for front-back-pairs
+    void print() const {
+    	std::cout << Form("TDCs\n trk:%1.4f\ntrkp:%1.4f\nq:%1.4f\ns800e1:%1.4f\ns800trg:%1.4f\nrf:%1.4f\nrfuw:%1.4f\ngt:%1.4f\nsi:%1.4f\nsiup:%1.4f\n-----\nevents_size:%lu\ntimestamp:%lld\nfound_de:%d, found_e:%d\n-------\n-------\n",
+    	tdc_trk, tdc_trkp, tdc_q, tdc_s800e1, tdc_s800trg, tdc_rf, tdc_rf_unwrap, tdc_gt, tdc_si, tdc_siup, events.size(), timestamp, found_de, found_e) << std::endl;
+    }
+};
+
+class trackingdet {
+public:
+    double timestamp=DEFAULT_NULL;
+    //TODO: Convert all to std::array
+    std::vector<float> xwires[MAXNWIRES_TRACK];
+    std::vector<float> ywires[MAXNWIRES_TRACK];
+    std::vector<float> xwiresf[MAXNWIRES_TRACK];
+    std::vector<float> ywiresf[MAXNWIRES_TRACK];
+    std::vector<float> xwiresf_nn[MAXNWIRES_TRACK];
+    std::vector<float> ywiresf_nn[MAXNWIRES_TRACK];
+
+    std::vector<float> xtimes[MAXNWIRES_TRACK];
+    std::vector<float> ytimes[MAXNWIRES_TRACK];
+    std::vector<float> xtimesf[MAXNWIRES_TRACK];
+    std::vector<float> ytimesf[MAXNWIRES_TRACK];
+    std::vector<float> cathode;
+
+    int multx; //how many x-wires fired?
+    int multy; //how many y-wires fired?
+    int multxf; //how many filtered x-wires fired?
+    int multyf; //how many filtered y-wires fired?
+
+    int multxt; //how many x-tdcwires fired?
+    int multyt; //how many y-tdcwires fired?
+
+    int multxtf; //how many x-tdcwires fired in window?
+    int multytf; //how many y-tdcwires fired in window?
+
+    float tot_cathode=0;
+    float tot_x=0;
+    float tot_y=0;
+    float tot_anode=0;
+
+    //list of x and y wires fired above energy threshold, within timing gate window
+    std::set<int> list_ywires;
+    std::set<int> list_xwires;
+
+    std::set<int> list_ytwires;
+    std::set<int> list_xtwires;
+
+    //position of the vertex estimated by up to 2 neighbouring wires firing together within window  
+    double xpos=DEFAULT_NULL;
+    double ypos=DEFAULT_NULL;
+    bool clean_event = false;
+    int maxnx=-124, maxny=-124;
+    int nnx = 0; //nearest neighbour x wires set this to +/- 1 if present w.r.t. wire maxnx
+    int nny = 0; //nearest neighbour y wires set this to +/- 1 if present w.r.t  wire maxny
+    bool clean_event_no_timing = false;
+    bool clean_single_xy_event = false;
+    bool clean_single_xy_event_no_timing = false;
+    void Reset() {
+        /***
+        Resets all data members.
+        **/
+        cathode.clear();
+        for(int i=0; i<MAXNWIRES_TRACK; i++) {
+            xwires[i].clear();
+            ywires[i].clear();
+            xwiresf[i].clear();
+            ywiresf[i].clear();
+            xwiresf_nn[i].clear();
+            ywiresf_nn[i].clear();
+            xtimes[i].clear();
+            ytimes[i].clear();
+            xtimesf[i].clear();
+            ytimesf[i].clear();
+        }
+        list_xwires.clear();
+        list_ywires.clear();
+        list_xtwires.clear();
+        list_ytwires.clear();
+
+        nnx = 0;
+        nny = 0;
+        clean_event_no_timing = false;
+        clean_event = false; //x, y both fire above thresh, xt, yt present within broad coinc window
+        xpos=DEFAULT_NULL;
+        ypos=DEFAULT_NULL;
+        maxnx=-124;
+        maxny=-124;
+        multx=0; //how many x-wires fired?
+        multy=0; //how many y-wires fired?
+        multxf=0; //how many filt x-wires fired?
+        multyf=0; //how many filt y-wires fired?
+        multxt=0; //how many tdc x-wires fired?
+        multyt=0; //how many tdc y-wires fired?
+        multxtf=0; //how many tdc x-wires fired?
+        multytf=0; //how many tdc y-wires fired?
+    }
+
+    trackingdet() {
+        Reset();
+    };
+    trackingdet(const std::vector<type19Raw>& orvec);
+};
+
+const float alpha = 0.0;
+int matchchantype(unsigned short chan, const std::array<orruba_params,MAX_ORRUBA_CHANS>& index, const std::string& label);
+void initialize_orruba(const std::string & filename1, const std::string& filename2, const std::string& filename3, const std::string& filename4);//,
+int parse_orruba_data(const unsigned short* buffer, int32_t length, type19Raw& oraw_event);
+#endif

gnuplot_lookerup

@@ -0,0 +1,80 @@
+
+# --- configuration ---
+datafilec = "../ANASEN_analysis/PW_study_Sudarsan/cathode1_2.dat"
+datafilea = "../ANASEN_analysis/PW_study_Sudarsan/anode1_2.dat"
+datafileqr = "../ANASEN_analysis/PW_study_Sudarsan/qqq_rings_geometry.dat"
+datafileqw = "../ANASEN_analysis/PW_study_Sudarsan/qqq_wedges_geometry.dat"
+
+set datafile separator ","
+
+N = 24   #number of cids
+
+# --- allocate arrays ---
+array ca1x[N]
+array an1x[N]
+array ca1y[N]
+array an1y[N]
+array ca1z[N]
+array an1z[N]
+array ca2x[N]
+array an2x[N]
+array ca2y[N]
+array an2y[N]
+array ca2z[N]
+array an2z[N]
+array qqqW_Phi[64] 
+array qqqR_Rho[16]
+
+# --- initialize (optional but safe) ---
+do for [i=1:N] {
+    ca1x[i] = ca1y[i] = ca1z[i] = NaN;
+    ca2x[i] = ca2y[i] = ca2z[i] = NaN;
+    an1x[i] = an1y[i] = an1z[i] = NaN;
+    an2x[i] = an2y[i] = an2z[i] = NaN;
+}
+
+do for [j=1:16] {
+    do for [i=1:4] {
+        qqqW_Phi[i*4+j] = NaN
+    }
+    qqqR_Rho[j] = NaN
+}
+
+
+# --- load table into arrays ---
+stats datafilec using \
+    (ca1x[int($1)+1] = $2, \
+     ca1y[int($1)+1] = $3, \
+     ca1z[int($1)+1] = $4, \
+     ca2x[int($1)+1] = $5, \
+     ca2y[int($1)+1] = $6, \
+     ca2z[int($1)+1] = $7, 0) \
+    nooutput
+
+stats datafilea using \
+    (an1x[int($1)+1] = $2, \
+     an1y[int($1)+1] = $3, \
+     an1z[int($1)+1] = $4, \
+     an2x[int($1)+1] = $5, \
+     an2y[int($1)+1] = $6, \
+     an2z[int($1)+1] = $7, 0) \
+    nooutput
+
+set datafile separator "\t"
+stats datafileqr using \
+    (qqqR_Rho[int($1)+1] = $2, 0) \
+    nooutput
+
+stats datafileqw using \
+    (qqqW_Phi[(int($1))*4+int($2)+1] = $3, 0) \
+    nooutput
+
+#do for [i=1:N] {
+#    print sprintf("Data[%d] = %g", i, an1x[i])
+#    print sprintf("Data[%d] = %g", i, ca1x[i])
+#}
+
+set datafile separator whitespace
+#plot '< cat /tmp/coords | grep q' u (ca1x[int($2)+1]):(ca1y[int($2)+1]):(ca1x[int($2)+1]):(ca1y[int($2)+1]) w vector
+
+#pause mouse key

live_plotter

@@ -0,0 +1,41 @@
+#columns:
+#1 is anodeindex, 2 is cathodeindex, 3,4,5 is crossover x,y,z, 6 is alpha - 7,8,9 are anode loci (first plane), 10, 11, 12 are cathode loci (first plane). 
+#13,14,15 are anode second plane loci, 16,17,18 are cathode second plane loci
+
+#set term qt font 'Verdana,10' size 1920,1080
+set term x11 font 'terminus-16' size 1920,1080 noraise
+#set key outside
+pause_status=0
+load 'gnuplot_lookerup'
+bind all "x" "pause_status=!pause_status"
+set datafile separator ","
+
+#anode-cathode-raw-loci, junction points colored same as anode wires
+set title sprintf('press ctrl+q to redraw, Ctrl+C the terminal to exit. pause\_status=%d',pause_status)
+set xlabel 'x-axis'
+set ylabel 'y-axis'
+#set view 212,110,1,1
+
+set xrange [-100:100]
+set yrange [-100:100]
+set zrange [-200:200]
+
+set zlabel 'z-axis'
+splot '../ANASEN_analysis/PW_study_Sudarsan/results_zsort.dat' u 7:8:9:($13-$7):($14-$8):($15-$9):1 with vectors nohead dt 1 lc variable title 'anodes',\
+	'' u 10:11:12:($16-$10):($17-$11):($18-$12):2 w vectors nohead dt 7 lc variable title 'cathodes',\
+    '../ANASEN_analysis/PW_study_Sudarsan/results_zsort.dat' u 3:4:5:1 w lp ls 7 lc variable title 'vertices colored anode',\
+    '< cat /tmp/coords | grep an' u (an1x[int($2)+1]):(an1y[int($2)+1]):(an1z[int($2)+1]):(an2x[int($2)+1]-an1x[int($2)+1]):(an2y[int($2)+1]-an1y[int($2)+1]):(an2z[int($2)+1]-an1z[int($2)+1]) w vector nohead lc 'black' lw 8 title 'selected\_anode',\
+    '< cat /tmp/coords | grep ca' u (ca1x[int($2)+1]):(ca1y[int($2)+1]):(ca1z[int($2)+1]):(ca2x[int($2)+1]-ca1x[int($2)+1]):(ca2y[int($2)+1]-ca1y[int($2)+1]):(ca2z[int($2)+1]-ca1z[int($2)+1]) w vectors nohead lc 'grey' dt 2 lw 8 title 'selected\_cath',\
+    '< cat /tmp/coords | grep q' u (qqqR_Rho[int($3)+1]*cos((pi/180.)*qqqW_Phi[int($2)*4+int($5)+1])):(qqqR_Rho[int($3)+1]*sin((pi/180.)*qqqW_Phi[int($2)*4+int($5)+1])):(128) w p ps 4 lc 'black' title 'qqqvertex'
+	#'' u 7:8:9:($1):1 w labels font ',10' tc variable notitle,\
+	#'' u 13:14:15:($1):1 w labels font ',10' tc variable notitle,\
+	#'' u 10:11:12:($2):2 w labels font ',10' tc variable notitle,\
+	#'' u 16:17:18:($2):2 w labels font ',10' tc variable notitle
+if(pause_status==0) {
+    pause 1.0
+} else {
+    reread
+    pause mouse close
+    pause_status=!pause_status
+}
+reread

pccal/anode_gainmatch.C

@@ -0,0 +1,67 @@
+void anode_gainmatch(){
+    TFile *f = new TFile("../results_run16.root");
+
+    TH2F *pc_index_h2d = (TH2F*)(f->Get("hRawPC/PC_Index_Vs_Energy"));
+    std::cout << pc_index_h2d << std::endl;
+    TCanvas c("c1","c1",0,0,1600,800);
+    //TCanvas c_g("cg","cg",0,900,400,400);
+    c.Divide(2,1);
+    auto c1=c.cd(1);
+    pc_index_h2d->Draw("COLZ");
+    pc_index_h2d->GetYaxis()->SetRangeUser(240,5000);
+    auto c2=c.cd(2);
+    c2->SetLogy();
+    TH1F *h_1d=NULL;
+    int bin_index=1;
+    std::vector<std::vector<double>> all_peaks;
+    std::vector<int> found_wire_list;
+    while(bin_index<=24) {
+        h_1d=(TH1F*)(pc_index_h2d->ProjectionY("_py",bin_index,bin_index));
+        auto c1 = c.cd(1);
+        TBox box(pc_index_h2d->GetXaxis()->GetBinLowEdge(bin_index),0,pc_index_h2d->GetXaxis()->GetBinUpEdge(bin_index),pc_index_h2d->GetYaxis()->GetXmax());
+        box.SetFillColorAlpha(kYellow+3,0.3);
+        box.Draw("SAME");
+        c1->Modified(); c1->Update();
+        //while(c1->WaitPrimitive());
+
+        TSpectrum s;
+        auto c2 = c.cd(2);
+        h_1d->Draw();
+        c2->Modified(); c2->Update();
+        int npeaks = s.Search(h_1d,8,"",0.02); std::cout << npeaks << std::endl;
+        if(npeaks>=3) {
+            std::vector<double> xpeaks(s.GetPositionX(),s.GetPositionX()+npeaks);
+            std::sort(xpeaks.begin(),xpeaks.end(),std::greater());
+            found_wire_list.push_back((int)pc_index_h2d->GetXaxis()->GetBinCenter(bin_index));
+            all_peaks.push_back(xpeaks);
+        }
+        while(c2->WaitPrimitive());
+        bin_index++;
+    }
+    c.cd(2)->SetLogy(kFALSE);
+    gStyle->SetOptFit(1111);
+
+    std::ofstream outfile("anode_gm_coeffs.dat");
+    outfile <<  found_wire_list.at(0) << " "
+            << 1.0  << " "
+            << 0.0  << std::endl;
+
+    for(int i=0; i<all_peaks.size(); i++){
+        if(i==1) continue;
+        TGraph g(all_peaks.at(i).size(), all_peaks.at(i).data(), all_peaks.at(1).data());
+        auto c2 = c.cd(2);
+        g.SetMarkerStyle(20);
+        //g.Print();
+        g.Draw("AP");
+        g.Fit("pol1");
+        outfile <<  found_wire_list.at(i) << " "
+                << ((TF1*)g.FindObject("pol1"))->GetParameter(1) << " "
+                << ((TF1*)g.FindObject("pol1"))->GetParameter(0) << std::endl;
+        c2->Modified();
+        c2->Update();
+        while(c2->WaitPrimitive());
+    }
+    outfile.close();
+    f->Close();
+    return;
+}

pccal/anode_gm_coeffs.dat

@@ -0,0 +1,24 @@
+0 1 0
+0 0.937314 -16.871
+2 0.965461 -1.54376
+3 0.926501 -3.27662
+4 0.905634 2.54577
+5 0.905634 -11.0387
+6 0.853919 6.23079
+7 0.945588 -9.54044
+8 0.884454 -11.8262
+9 0.922501 -3.42538
+10 0.903053 9.28069
+11 0.914653 9.87642
+12 0.965332 13.2526
+13 0.923847 -3.41775
+14 0.93845 25.9901
+15 0.955424 12.324
+16 0.95116 4.99595
+17 0.910745 2.86648
+18 0.941376 4.57217
+19 0.871622 932.111
+20 1.00624 7.86358
+21 0.969834 -45.001
+22 0.89304 -31.5635
+23 0.933226 4.02193

pccal/cathode_gainmatch.C

@@ -0,0 +1,69 @@
+void cathode_gainmatch(){
+    TFile *f = new TFile("../results_run17.root");
+    TH2F *pc_index_h2d = (TH2F*)(f->Get("hRawPC/PC_Index_Vs_Energy"));
+    std::cout << pc_index_h2d << std::endl;
+    TCanvas c("c1","c1",0,0,1600,800);
+    //TCanvas c_g("cg","cg",0,900,400,400);
+    c.Divide(2,1);
+    auto c1=c.cd(1);
+    pc_index_h2d->Draw("COLZ");
+    pc_index_h2d->GetYaxis()->SetRangeUser(600,pc_index_h2d->GetYaxis()->GetXmax());
+    auto c2=c.cd(2);
+    c2->SetLogy();
+    TH1F *h_1d=NULL;
+    int bin_index=25;
+    std::vector<double> pulser_heights = {0.01,0.05,0.1,0.15,0.2,0.25,0.3,0.5};
+    std::vector<std::vector<double>> all_peaks;
+    std::vector<int> found_wire_list;
+    while(bin_index<=48) {
+        h_1d=(TH1F*)(pc_index_h2d->ProjectionY("_py",bin_index,bin_index));
+        auto c1 = c.cd(1);
+        TBox box(pc_index_h2d->GetXaxis()->GetBinLowEdge(bin_index),0,pc_index_h2d->GetXaxis()->GetBinUpEdge(bin_index),pc_index_h2d->GetYaxis()->GetXmax());
+        box.SetFillColorAlpha(kYellow+3,0.3);
+        box.Draw("SAME");
+        c1->Modified(); c1->Update();
+        //while(c1->WaitPrimitive());
+
+        TSpectrum s;
+        auto c2 = c.cd(2);
+        h_1d->Draw();
+        c2->Modified(); c2->Update();
+        int npeaks = s.Search(h_1d,20,"",0.1); std::cout << npeaks << std::endl;
+        if(npeaks==8) {
+            std::vector<double> xpeaks(s.GetPositionX(),s.GetPositionX()+npeaks);
+            for(int i=0; i<8; i++) {
+                std::cout << pc_index_h2d->GetXaxis()->GetBinCenter(bin_index) << " " << xpeaks.at(i) << " " << xpeaks.at(i)/pulser_heights.at(i) << std::endl;
+            }
+            std::sort(xpeaks.begin(),xpeaks.end(),std::greater());
+            found_wire_list.push_back((int)pc_index_h2d->GetXaxis()->GetBinCenter(bin_index));
+            all_peaks.push_back(xpeaks);
+        }
+        while(c2->WaitPrimitive());
+        bin_index++;
+    }
+    c.cd(2)->SetLogy(kFALSE);
+    gStyle->SetOptFit(1111);
+
+    std::ofstream outfile("cathode_gm_coeffs.dat");
+    outfile <<  found_wire_list.at(0) << " "
+            << 1.0  << " "
+            << 0.0  << std::endl;
+
+    for(int i=1; i<all_peaks.size(); i++){
+        TGraph g(all_peaks.at(i).size(), all_peaks.at(i).data(), all_peaks.at(0).data());
+        auto c2 = c.cd(2);
+        g.SetMarkerStyle(20);
+        //g.Print();
+        g.Draw("AP");
+        g.Fit("pol1");
+        outfile <<  found_wire_list.at(i) << " "
+                << ((TF1*)g.FindObject("pol1"))->GetParameter(1) << " "
+                << ((TF1*)g.FindObject("pol1"))->GetParameter(0) << std::endl;
+        c2->Modified();
+        c2->Update();
+        while(c2->WaitPrimitive());
+    }
+    outfile.close();
+    f->Close();
+    return;
+}

pccal/cathode_gm_coeffs.dat

@@ -0,0 +1,21 @@
+24 1 0
+25 0.941896 6.16135
+26 0.980284 2.86886
+27 0.983166 -3.82952
+28 0.978704 -2.89713
+29 0.964947 2.25786
+30 0.94514 0.925074
+31 0.977231 1.6493
+32 0.919527 5.82742
+33 0.972243 2.88061
+34 0.928892 7.61384
+35 0.947376 -0.644223
+36 0.875342 6.066
+38 0.970953 6.262
+40 0.918408 -3.27891
+41 0.913619 4.11288
+42 0.954083 2.21261
+43 0.993037 5.48924
+45 0.926406 -19.719
+46 1.00459 5.14574
+47 0.942483 5.54183

pccal/pc_gm_coeffs.dat

@@ -0,0 +1,49 @@
+#Histogram Number	Slope	Intercept
+0 0.937314 -16.871
+1 0 0
+2 0.965461 -1.54376
+3 0.926501 -3.27662
+4 0.905634 2.54577
+5 0.905634 -11.0387
+6 0.853919 6.23079
+7 0.945588 -9.54044
+8 0.884454 -11.8262
+9 0.922501 -3.42538
+10 0.903053 9.28069
+11 0.914653 9.87642
+12 0.965332 13.2526
+13 0.923847 -3.41775
+14 0.93845 25.9901
+15 0.955424 12.324
+16 0.95116 4.99595
+17 0.910745 2.86648
+18 0.941376 4.57217
+19 0.871622 932.111
+20 1.00624 7.86358
+21 0.969834 -45.001
+22 0.89304 -31.5635
+23 0.933226 4.02193
+24 0 0
+25 0.941896 6.16135
+26 0.980284 2.86886
+27 0.983166 -3.82952
+28 0.978704 -2.89713
+29 0.964947 2.25786
+30 0.94514 0.925074
+31 0.977231 1.6493
+32 0.919527 5.82742
+33 0.972243 2.88061
+34 0.928892 7.61384
+35 0.947376 -0.644223
+36 0.875342 6.066
+37 0    0
+38 0.970953 6.262
+39 0    0
+40 0.918408 -3.27891
+41 0.913619 4.11288
+42 0.954083 2.21261
+43 0.993037 5.48924
+44  0   0
+45 0.926406 -19.719
+46 1.00459 5.14574
+47 0.942483 5.54183

pccal/slope_intercept_26Al.dat

@@ -0,0 +1,50 @@
+#Histogram Number	Slope	Intercept
+#Histogram Number	Slope	Intercept
+0	0.931015	-1.35431
+1	1	-1.87356e-10
+2	0.964185	1.49989
+3	0.92638	-1.30621
+4	0.905569	1.00834
+5	0.901182	0.470903
+6	0.853932	3.32687
+7	0.942785	1.08887
+8	0.878904	-0.0107433
+9	0.922662	-2.32259
+10	0.903343	8.38332
+11	0.914227	6.56108
+12	0.961008	23.0982
+13	0.920976	5.22104
+14	0.936584	31.5073
+15	0.959044	5.43267
+16	0.95263	-0.404053
+17	0.90953	4.82833
+18	0.940277	10.3629
+19	0.86746	-17.8678
+20	1.00683	4.76371
+21	0.968342	-43.9496
+22	0.892882	-32.0742
+23	0.933615	1.10704
+24	1	-2.89219e-10
+25	0.942098	-0.105169
+26	0.980862	-0.732032
+27	0.982975	-2.22704
+28	0.978815	-1.51477
+29	0.965245	-2.19515
+30	0.945384	-0.892599
+31	0.977408	-0.908592
+32	0.919546	3.25464
+33	0.972194	2.44956
+34	0.92852	5.44745
+35	0.947098	1.40531
+36	0.875491	-1.13145
+37	1	0
+38	0.970862	2.86019
+39	1	0
+40	0.91793	-3.80615
+41	0.913897	-2.12964
+42	0.954014	-0.760604
+43	0.993616	-1.40278
+44	1	0
+45	0.926169	-21.2016
+46	1.00577	-2.14281
+47	0.943312	-1.26464

run_sx3.sh

@@ -7,10 +7,11 @@
 #root -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_004_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run04.root;
 
 #alpha+gas 27Al
-#root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run12.root;
+root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run09.root;
+root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run12.root;
 
 #protons+gas, 27Al
-#root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_022_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run22.root;
+root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_022_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run22.root;
 #root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_021_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run21.root;
 
 #27Al reaction data
@@ -20,11 +21,15 @@
 
 #root -l -x  results_run19.root results_run12.root -e "new TBrowser"
 
+#17F pulser runs
+#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/PulserRun_015_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run15.root;
+#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/PulserRun_016_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run16.root;
+#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/PulserRun_017_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run17.root;
+
 #17F alpha run with gas
-#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run18.root;
+root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_018_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run18.root;
 root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_019_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run19.root;
-#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run20.root;
+root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_020_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run20.root;
-#root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run21.root;
+root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run21.root;
 #root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Run_104_mapped.root -e 'tree->Process("MakeVertex.C+")'; mv Analyzer_SX3.root results_run104.root;
 #mv Analyzer_SX3.root results_run19.root;
-

scratch/25Mar26_17F_alpha+gaus/compare.C

@@ -0,0 +1,41 @@
+#include <signal.h>
+int quit=0;
+void handler(int signal){ quit=1;}
+void compare() {
+    signal(SIGINT,handler);
+    TCanvas c;
+    c.SetLogz();
+    TFile f18("../../results_run18.root");
+    TFile f19("../../results_run19.root");
+    TFile f20("../../results_run20.root");
+    TFile f21("../../results_run21.root");
+    TH2F *h18 = (TH2F*)(f18.Get("pcz_vs_sx3z"));
+    TH2F *h19 = (TH2F*)(f19.Get("pcz_vs_sx3z"));
+    TH2F *h20 = (TH2F*)(f20.Get("pcz_vs_sx3z"));
+    TH2F *h21 = (TH2F*)(f21.Get("pcz_vs_sx3z"));
+
+    h18->GetYaxis()->SetRangeUser(-200,200);
+    h19->GetYaxis()->SetRangeUser(-200,200);
+    h20->GetYaxis()->SetRangeUser(-200,200);
+    h21->GetYaxis()->SetRangeUser(-200,200);
+    h18->GetXaxis()->SetRangeUser(0,100);
+    h19->GetXaxis()->SetRangeUser(0,100);
+    h20->GetXaxis()->SetRangeUser(0,100);
+    h21->GetXaxis()->SetRangeUser(0,100);
+
+    while(!quit) {
+        h18->Draw("COLZ");
+        gPad->Modified(); gPad->Update(); while(gPad->WaitPrimitive());
+        h19->Draw("COLZ");
+        gPad->Modified(); gPad->Update(); while(gPad->WaitPrimitive());
+        h20->Draw("COLZ");
+        gPad->Modified(); gPad->Update(); while(gPad->WaitPrimitive());
+        h21->Draw("COLZ");
+        gPad->Modified(); gPad->Update(); while(gPad->WaitPrimitive());
+    }
+
+    f18.Close();
+    f19.Close();
+    f20.Close();
+    f21.Close();
+}

scratch/25Mar26_17F_alpha+gaus/fit.log

@@ -0,0 +1,279 @@
+
+
+*******************************************************************************
+Wed Mar 25 17:10:21 2026
+
+
+FIT:    data read from 'first_test.txt'
+        format = z
+        #datapoints = 4
+        residuals are weighted equally (unit weight)
+
+function used for fitting: A*x+B
+fitted parameters initialized with current variable values
+
+Warning: Initial value of parameter 'B' is zero.
+iter      chisq       delta/lim  lambda   A             B            
+   0 3.7754834810e+03   0.00e+00  2.97e+01    1.000000e+00   1.000000e-30
+BREAK: Singular matrix in Invert_RtR
+
+*******************************************************************************
+Wed Mar 25 17:10:29 2026
+
+
+FIT:    data read from 'first_test.txt'
+        format = z
+        #datapoints = 4
+        residuals are weighted equally (unit weight)
+
+function used for fitting: A*x+B
+fitted parameters initialized with current variable values
+
+iter      chisq       delta/lim  lambda   A             B            
+   0 3.5907054810e+03   0.00e+00  2.97e+01    1.000000e+00   1.000000e+00
+   5 9.9465658244e+01  -1.35e-03  2.97e-04    1.682744e+00   4.503012e+01
+
+After 5 iterations the fit converged.
+final sum of squares of residuals : 99.4657
+rel. change during last iteration : -1.34817e-08
+
+degrees of freedom    (FIT_NDF)                        : 2
+rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 7.05215
+variance of residuals (reduced chisquare) = WSSR/ndf   : 49.7328
+
+Final set of parameters            Asymptotic Standard Error
+=======================            ==========================
+A               = 1.68274          +/- 0.1265       (7.517%)
+B               = 45.0301          +/- 5.311        (11.79%)
+
+correlation matrix of the fit parameters:
+                A      B      
+A               1.000 
+B               0.748  1.000 
+
+
+*******************************************************************************
+Wed Mar 25 17:14:06 2026
+
+
+FIT:    data read from 'first_test.txt'
+        format = z
+        #datapoints = 4
+        residuals are weighted equally (unit weight)
+
+function used for fitting: A*x+B
+fitted parameters initialized with current variable values
+
+Warning: Initial value of parameter 'B' is zero.
+iter      chisq       delta/lim  lambda   A             B            
+   0 3.7754834810e+03   0.00e+00  3.38e+01    1.000000e+00   1.000000e-30
+BREAK: Singular matrix in Invert_RtR
+
+*******************************************************************************
+Wed Mar 25 17:14:13 2026
+
+
+FIT:    data read from 'first_test.txt'
+        format = z
+        #datapoints = 4
+        residuals are weighted equally (unit weight)
+
+function used for fitting: A*x+B
+fitted parameters initialized with current variable values
+
+iter      chisq       delta/lim  lambda   A             B            
+   0 3.9682614810e+03   0.00e+00  3.38e+01    1.000000e+00   1.000000e+00
+   5 3.4734096263e+01  -3.36e-04  3.38e-04    5.876258e-01  -2.681171e+01
+
+After 5 iterations the fit converged.
+final sum of squares of residuals : 34.7341
+rel. change during last iteration : -3.36116e-09
+
+degrees of freedom    (FIT_NDF)                        : 2
+rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 4.16738
+variance of residuals (reduced chisquare) = WSSR/ndf   : 17.367
+
+Final set of parameters            Asymptotic Standard Error
+=======================            ==========================
+A               = 0.587626         +/- 0.04417      (7.517%)
+B               = -26.8117         +/- 2.112        (7.877%)
+
+correlation matrix of the fit parameters:
+                A      B      
+A               1.000 
+B               0.163  1.000 
+
+
+*******************************************************************************
+Wed Mar 25 17:14:34 2026
+
+
+FIT:    data read from 'first_test.txt'
+        format = z
+        #datapoints = 4
+        residuals are weighted equally (unit weight)
+
+function used for fitting: A*x+B
+fitted parameters initialized with current variable values
+
+iter      chisq       delta/lim  lambda   A             B            
+   0 3.9682614810e+03   0.00e+00  3.38e+01    1.000000e+00   1.000000e+00
+   5 3.4734096263e+01  -3.36e-04  3.38e-04    5.876258e-01  -2.681171e+01
+
+After 5 iterations the fit converged.
+final sum of squares of residuals : 34.7341
+rel. change during last iteration : -3.36116e-09
+
+degrees of freedom    (FIT_NDF)                        : 2
+rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 4.16738
+variance of residuals (reduced chisquare) = WSSR/ndf   : 17.367
+
+Final set of parameters            Asymptotic Standard Error
+=======================            ==========================
+A               = 0.587626         +/- 0.04417      (7.517%)
+B               = -26.8117         +/- 2.112        (7.877%)
+
+correlation matrix of the fit parameters:
+                A      B      
+A               1.000 
+B               0.163  1.000 
+
+
+*******************************************************************************
+Wed Mar 25 18:00:20 2026
+
+
+FIT:    data read from 'first_test.txt'
+        format = z
+        #datapoints = 5
+        residuals are weighted equally (unit weight)
+
+function used for fitting: A*x+B
+fitted parameters initialized with current variable values
+
+iter      chisq       delta/lim  lambda   A             B            
+   0 4.3038838810e+03   0.00e+00  3.03e+01    1.000000e+00   1.000000e+00
+   5 7.7126310821e+01  -2.27e-05  3.03e-04    5.892187e-01  -2.534379e+01
+
+After 5 iterations the fit converged.
+final sum of squares of residuals : 77.1263
+rel. change during last iteration : -2.27353e-10
+
+degrees of freedom    (FIT_NDF)                        : 3
+rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 5.07038
+variance of residuals (reduced chisquare) = WSSR/ndf   : 25.7088
+
+Final set of parameters            Asymptotic Standard Error
+=======================            ==========================
+A               = 0.589219         +/- 0.05373      (9.119%)
+B               = -25.3438         +/- 2.301        (9.08%)
+
+correlation matrix of the fit parameters:
+                A      B      
+A               1.000 
+B               0.171  1.000 
+
+
+*******************************************************************************
+Wed Mar 25 18:00:37 2026
+
+
+FIT:    data read from 'first_test.txt'
+        format = z
+        #datapoints = 4
+        residuals are weighted equally (unit weight)
+
+function used for fitting: A*x+B
+fitted parameters initialized with current variable values
+
+iter      chisq       delta/lim  lambda   A             B            
+   0 3.9682614810e+03   0.00e+00  3.38e+01    1.000000e+00   1.000000e+00
+   5 3.4734096263e+01  -3.36e-04  3.38e-04    5.876258e-01  -2.681171e+01
+
+After 5 iterations the fit converged.
+final sum of squares of residuals : 34.7341
+rel. change during last iteration : -3.36116e-09
+
+degrees of freedom    (FIT_NDF)                        : 2
+rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 4.16738
+variance of residuals (reduced chisquare) = WSSR/ndf   : 17.367
+
+Final set of parameters            Asymptotic Standard Error
+=======================            ==========================
+A               = 0.587626         +/- 0.04417      (7.517%)
+B               = -26.8117         +/- 2.112        (7.877%)
+
+correlation matrix of the fit parameters:
+                A      B      
+A               1.000 
+B               0.163  1.000 
+
+
+*******************************************************************************
+Wed Mar 25 18:00:50 2026
+
+
+FIT:    data read from 'first_test.txt'
+        format = z
+        #datapoints = 5
+        residuals are weighted equally (unit weight)
+
+function used for fitting: A*x+B
+fitted parameters initialized with current variable values
+
+iter      chisq       delta/lim  lambda   A             B            
+   0 4.3038838810e+03   0.00e+00  3.03e+01    1.000000e+00   1.000000e+00
+   5 7.7126310821e+01  -2.27e-05  3.03e-04    5.892187e-01  -2.534379e+01
+
+After 5 iterations the fit converged.
+final sum of squares of residuals : 77.1263
+rel. change during last iteration : -2.27353e-10
+
+degrees of freedom    (FIT_NDF)                        : 3
+rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 5.07038
+variance of residuals (reduced chisquare) = WSSR/ndf   : 25.7088
+
+Final set of parameters            Asymptotic Standard Error
+=======================            ==========================
+A               = 0.589219         +/- 0.05373      (9.119%)
+B               = -25.3438         +/- 2.301        (9.08%)
+
+correlation matrix of the fit parameters:
+                A      B      
+A               1.000 
+B               0.171  1.000 
+
+
+*******************************************************************************
+Wed Mar 25 18:01:05 2026
+
+
+FIT:    data read from 'first_test.txt'
+        format = z
+        #datapoints = 5
+        residuals are weighted equally (unit weight)
+
+function used for fitting: A*x+B
+fitted parameters initialized with current variable values
+
+iter      chisq       delta/lim  lambda   A             B            
+   0 4.3038838810e+03   0.00e+00  3.03e+01    1.000000e+00   1.000000e+00
+   5 7.7126310821e+01  -2.27e-05  3.03e-04    5.892187e-01  -2.534379e+01
+
+After 5 iterations the fit converged.
+final sum of squares of residuals : 77.1263
+rel. change during last iteration : -2.27353e-10
+
+degrees of freedom    (FIT_NDF)                        : 3
+rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 5.07038
+variance of residuals (reduced chisquare) = WSSR/ndf   : 25.7088
+
+Final set of parameters            Asymptotic Standard Error
+=======================            ==========================
+A               = 0.589219         +/- 0.05373      (9.119%)
+B               = -25.3438         +/- 2.301        (9.08%)
+
+correlation matrix of the fit parameters:
+                A      B      
+A               1.000 
+B               0.171  1.000 

scratch/25Mar26_17F_alpha+gaus/gnuplotter_first

@@ -0,0 +1,7 @@
+set term qt enhanced
+set key top left
+A=1;
+B=1;
+fit A*x+B 'first_test.txt' via A,B
+plot 'first_test.txt' w p, A*x+B w l title 'fit'
+pause mouse close

scratch/lise.dat

@@ -0,0 +1,16 @@
+#ein 1cm-eloss 3cm-eloss
+6   0.1446   0.4458
+5.5 0.1579   0.485
+5.  0.1657   0.5149
+4.5 0.1829   0.565
+4.  0.1976   0.617
+3.5 0.2184   0.6878
+3.  0.2441   0.780
+2.5 0.2774   0.911
+2.  0.3231   1.086
+1.5 0.3818   1.19
+1.  0.4242   0.9915
+0.5 0.3439   0.5
+0.2 0.1934   0.2
+0.1 0.1      0.1
+

scratch/phiscan.sh

@@ -0,0 +1,11 @@
+declare -i i=-7
+while [[ $i -lt 3 ]]; do
+    echo "offset_c1 = -4*k  +"$i"*k - TMath::TwoPi()/48;" > testing.h
+    cd ..
+    bash run_sx3.sh
+    j=$(echo $i+3 | bc)
+    cp results_run12.root results_run12case$j.root
+    cd -
+    i=i+1
+done
+echo "" > testing.h

scratch/scan_phisearch.C

@@ -0,0 +1,20 @@
+{
+    for(int i=-4; i<=5; i++) {
+        TFile *f = new TFile(Form("../results_run12case%d.root",i));
+        TH2F *h2 = (TH2F*)(f->Get("phicut/pczguess_vs_pc_int"));
+//        TH2F *h2 = (TH2F*)(f->Get("pcz_vs_sx3pczguess"));
+//        TH2F *h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ2"));
+        h2->SetTitle(Form("case%d",i));
+        h2->Draw("colz");
+        gPad->Modified();
+        gPad->Update();
+        while(gPad->WaitPrimitive());
+
+        h2 = (TH2F*)(f->Get("hPCQQQ/PC_XY_Projection_QQQ3"));
+        h2->Draw("colz");
+        gPad->Modified();
+        gPad->Update();
+        while(gPad->WaitPrimitive());
+        f->Close();
+    }
+}

scratch/test.dat

@@ -0,0 +1,7 @@
+1.00490190125271	1.81526549142349
+2.1078430941584	1.4359987348238
+3.02287578753202	1.14680783291655
+3.87254900280751	0.971396957989193
+4.67320260950942	0.852876096551793
+5.44117647716228	0.819690255349321
+6.02124184528305	0.800726917519336

scratch/testing.h

@@ -0,0 +1 @@
+

scratch/testme.dat

@@ -0,0 +1,29 @@
+{
+//========= Macro generated from object: Graph/Graph
+//========= by ROOT version6.36.04
+   
+   std::vector<Double_t> graph_x_vect0{ 1.00490190125271, 2.107843094158404, 3.022875787532016, 3.872549002807514, 4.673202609509424, 5.441176477162277, 6.02124184528305 };
+   std::vector<Double_t> graph_y_vect1{ 18152.65491423486, 14359.98734823804, 11468.07832916546, 9713.969579891931, 8528.760965517926, 8196.902553493204, 8007.269175193363 };
+   TGraph *graph = new TGraph(7, graph_x_vect0.data(), graph_y_vect1.data());
+   graph->SetName("Graph");
+   graph->SetTitle("Graph");
+   graph->SetFillStyle(1000);
+   
+   TH1F *Graph_histogram1 = new TH1F("Graph_histogram1", "Graph", 100, 0.4542482982760897, 6.473856231112497);
+   Graph_histogram1->SetMinimum(7324.589013313936);
+   Graph_histogram1->SetMaximum(19499.05190016373);
+   Graph_histogram1->SetDirectory(nullptr);
+   Graph_histogram1->SetStats(0);
+   Graph_histogram1->SetLineColor(TColor::GetColor("#000099"));
+   Graph_histogram1->GetXaxis()->SetLabelFont(42);
+   Graph_histogram1->GetXaxis()->SetTitleOffset(1);
+   Graph_histogram1->GetXaxis()->SetTitleFont(42);
+   Graph_histogram1->GetYaxis()->SetLabelFont(42);
+   Graph_histogram1->GetYaxis()->SetTitleFont(42);
+   Graph_histogram1->GetZaxis()->SetLabelFont(42);
+   Graph_histogram1->GetZaxis()->SetTitleOffset(1);
+   Graph_histogram1->GetZaxis()->SetTitleFont(42);
+   graph->SetHistogram(Graph_histogram1);
+   
+   graph->Draw();
+}

slope_intercept_results.dat

@@ -1,4 +1,4 @@
-Histogram Number	Slope	Intercept
+#Histogram Number	Slope	Intercept
 0	0.931015	-1.35431
 1	1	-1.87356e-10
 2	0.964185	1.49989

slope_intercept_results_17F.dat

@@ -0,0 +1,49 @@
+#Histogram Number	Slope	Intercept
+0 0.937314 -16.871
+1 0 0
+2 0.965461 -1.54376
+3 0.926501 -3.27662
+4 0.905634 2.54577
+5 0.905634 -11.0387
+6 0.853919 6.23079
+7 0.945588 -9.54044
+8 0.884454 -11.8262
+9 0.922501 -3.42538
+10 0.903053 9.28069
+11 0.914653 9.87642
+12 0.965332 13.2526
+13 0.923847 -3.41775
+14 0.93845 25.9901
+15 0.955424 12.324
+16 0.95116 4.99595
+17 0.910745 2.86648
+18 0.941376 4.57217
+19 0.871622 932.111
+20 1.00624 7.86358
+21 0.969834 -45.001
+22 0.89304 -31.5635
+23 0.933226 4.02193
+24 0 0
+25 0.941896 6.16135
+26 0.980284 2.86886
+27 0.983166 -3.82952
+28 0.978704 -2.89713
+29 0.964947 2.25786
+30 0.94514 0.925074
+31 0.977231 1.6493
+32 0.919527 5.82742
+33 0.972243 2.88061
+34 0.928892 7.61384
+35 0.947376 -0.644223
+36 0.875342 6.066
+37 0    0
+38 0.970953 6.262
+39 0    0
+40 0.918408 -3.27891
+41 0.913619 4.11288
+42 0.954083 2.21261
+43 0.993037 5.48924
+44  0   0
+45 0.926406 -19.719
+46 1.00459 5.14574
+47 0.942483 5.54183