modified: .vscode/c_cpp_properties.json

.vscode/c_cpp_properties.json

@@ -1,23 +1,11 @@
 {
     "configurations": [
         {
-            "name": "splitpole",
+            "name": "Hades",
             "includePath": [
                 "${workspaceFolder}/**",
-                "/home/splitpole/cern/root/**"
+                "/usr/include/x86_64-linux-gnu/qt6/**",
-            ],
+                "/usr/local/cern/root_v6.26.06/include/**"
-            "defines": [],
-            "compilerPath": "/usr/bin/gcc",
-            "cStandard": "c17",
-            "cppStandard": "gnu++17",
-            "intelliSenseMode": "linux-gcc-x64"
-        },
-        {
-            "name": "Ryan",
-            "includePath": [
-                "${workspaceFolder}/**",
-                "/home/ryan/Downloads/root_build/**",
-                "/home/ryan/Downloads/root_build/include"
             ],
             "defines": [],
             "compilerPath": "/usr/bin/gcc",
@@ -29,7 +17,8 @@
             "name": "RyanUbuntu",
             "includePath": [
                 "${workspaceFolder}/**",
-                "/opt/root/**"
+                "/usr/include/x86_64-linux-gnu/qt6/**",
+                "/opt/root/include/**"
             ],
             "defines": [],
             "compilerPath": "/usr/bin/gcc",
@@ -38,10 +27,11 @@
             "intelliSenseMode": "linux-gcc-x64"
         },
         {
-            "name": "RyanHome",
+            "name": "Anasen",
             "includePath": [
                 "${workspaceFolder}/**",
-                "/home/ryan/root_v6.30.06/**"
+                "/usr/include/x86_64-linux-gnu/qt6/**",
+                "/opt/root/include/**"
             ],
             "defines": [],
             "compilerPath": "/usr/bin/gcc",
@@ -50,10 +40,26 @@
             "intelliSenseMode": "linux-gcc-x64"
         },
         {
-            "name": "Dirac",
+            "name": "Splitpole",
             "includePath": [
                 "${workspaceFolder}/**",
-                "/usr/opt/root/**"
+                "/usr/include/x86_64-linux-gnu/qt6/**",
+                "/home/splitpole/cern/root/include/**",
+                "/usr/include/x86_64-linux-gnu/qt6/QtWidgets",
+                "/usr/include/x86_64-linux-gnu/qt6/QtCore"
+            ],
+            "defines": [],
+            "compilerPath": "/usr/bin/gcc",
+            "cStandard": "c17",
+            "cppStandard": "gnu++17",
+            "intelliSenseMode": "linux-gcc-x64"
+        },
+        {
+            "name": "Penguin",
+            "includePath": [
+                "${workspaceFolder}/**",
+                "/usr/include/x86_64-linux-gnu/qt6/**",
+                "/usr/local/cern/root/include/**",
             ],
             "defines": [],
             "compilerPath": "/usr/bin/gcc",

Analysis.C

@@ -16,5 +16,5 @@ void Analysis(int start, int end) {
 
 // Define a macro with the same name as the script
 void Analysis() {
-    Analysis(72, 194); // Adjust the range if needed
+    Analysis(150, 194); // Adjust the range if needed
 }

Analyzer.C

@@ -58,18 +58,17 @@ void Analyzer::Begin(TTree * /*tree*/){
   hqqqVpcIndex->SetNdivisions( -612, "x");
   hqqqVpcIndex->SetNdivisions( -12, "y");
 
-  hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 400, 0, 5000);
+  hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 800, 0, 16000);
   hqqqVpcE->SetNdivisions( -612, "x");
   hqqqVpcE->SetNdivisions( -12, "y");
   
-  hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000);
+  hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 800, 0, 16000);
   hsx3VpcE->SetNdivisions( -612, "x");
   hsx3VpcE->SetNdivisions( -12, "y");
 
   hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
 
   hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000);
-
   sx3_contr.ConstructGeo();
   pw_contr.ConstructGeo();
 
@@ -283,7 +282,6 @@ Bool_t Analyzer::Process(Long64_t entry){
 
     float aE = 0;
     float cE = 0;
-    bool multi_an =false;
     // if( ID[0].first < 1 ) {
     //   aID = pc.ch[ID[0].second];
     //   cID = pc.ch[ID[1].second];
@@ -293,21 +291,28 @@ Bool_t Analyzer::Process(Long64_t entry){
     // }
     // printf("anode= %d, cathode = %d\n", aID, cID);
 
-  // for( int k = 0; k < qqq.multi; k++){
+   for( int k = 0; k < qqq.multi; k++){
-  //   if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){
+    if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){
+    
+      int multi_an =0;
       for(int l=0;l<E.size();l++){
-        if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){
+      if(E[l].first<24){
-          if(!multi_an){
+          multi_an++;
+        }
+      }
+
+      if(multi_an>=1){
+        for(int l=0;l<E.size();l++){
+          if(E[l].first<24 && E[l].first!=19 && E[l].first!=12){
             aE = E[l].second;
-          }
+          }else if(E[l].first>24){
-          multi_an=true;
+            cE = E[l].second;
-        }
+          }
-        else {
+        }
-          cE = E[l].second + cE;
+      }
+    
      }
    }
-    // }
-  // }
     hanVScatsum->Fill(aE,cE);
       
     if( ID[0].first < 1 ) {

Armory/#ClassPW.h#

@@ -0,0 +1,283 @@
+#ifndef ClassPW_h
+#define ClassPW_h
+
+#include <cstdio>
+#include <TMath.h>
+#include <TVector3.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;
+  }
+};
+
+//!########################################################
+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;}
+
+  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, PWHitInfo 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 | 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
+    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 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 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose){
+
+  trackPos = sx3Pos;
+
+  double p1 = TMath::Abs(hitInfo.nearestDist.first + gRandom->Gaus(0, sigmaA));
+  double p2 = TMath::Abs(hitInfo.nextNearestDist.first + gRandom->Gaus(0, sigmaA));
+  double fracA = p1 / (p1 + p2);
+  short anodeID1 = hitInfo.nearestWire.first;
+  short anodeID2 = hitInfo.nextNearestWire.first;
+  TVector3 shiftA1 = (An[anodeID2].first - An[anodeID1].first) * fracA;
+  TVector3 shiftA2 = (An[anodeID2].second - An[anodeID1].second) * fracA;
+
+  double q1 = TMath::Abs(hitInfo.nearestDist.second + gRandom->Gaus(0, sigmaC));
+  double q2 = TMath::Abs(hitInfo.nextNearestDist.second + gRandom->Gaus(0, sigmaC));
+  double fracC = q1 / (q1 + q2);
+  short cathodeID1 = hitInfo.nearestWire.second;
+  short cathodeID2 = hitInfo.nextNearestWire.second;
+  TVector3 shiftC1 = (Ca[cathodeID2].first - Ca[cathodeID1].first) * fracC;
+  TVector3 shiftC2 = (Ca[cathodeID2].second - Ca[cathodeID1].second) * fracC;
+
+  TVector3 a1 = An[anodeID1].first + shiftA1;
+  TVector3 a2 = An[anodeID1].second + shiftA2;
+
+  TVector3 c1 = Ca[cathodeID1].first + shiftC1;
+  TVector3 c2 = Ca[cathodeID1].second + shiftC2;
+
+  TVector3 n1 = (a1 - a2).Cross((sx3Pos - a2)).Unit();
+  TVector3 n2 = (c1 - c2).Cross((sx3Pos - c2)).Unit();
+
+  // if the handiness of anode and cathode revered, it should be n2 cross n1
+  trackVec = (n2.Cross(n1)).Unit();
+
+  if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg()); 
+
+}
+
+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 trackPos.Z() - rho / TMath::Tan(theta);
+
+}
+
+#endif 

Armory/ClassPW.h

@@ -4,6 +4,7 @@
 #include <cstdio>
 #include <TMath.h>
 #include <TVector3.h>
+#include <TRandom.h>
 
 struct PWHitInfo{
   std::pair<short, short>   nearestWire; // anode, cathode
@@ -65,8 +66,6 @@ public:
   void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
   void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false);
 
-  double CircularMean(std::vector<std::pair<int, double>> wireList);
-
   void Print(){
     printf("     The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first, 
                                                                          hitInfo.nearestDist.first, 
@@ -282,21 +281,4 @@ inline double PW::GetZ0(){
 
 }
 
-inline double PW::CircularMean(std::vector<std::pair<int, double>> wireList){
-
-  //use unit vector, wireID start from Zero
-  double xCom = 0, yCom = 0;
-  for( size_t i = 0; i < wireList.size() ; i++){
-    xCom += TMath::Cos(TMath::TwoPi() * wireList[i].first / nWire) * wireList[i].second;
-    yCom += TMath::Sin(TMath::TwoPi() * wireList[i].first / nWire) * wireList[i].second;
-  }
-
-  //calculate the angle of the summed unit vectors
-  double angle = TMath::ATan2(yCom, xCom);
-  if( angle < 0 ) angle += TMath::TwoPi(); // convert the angle from 0 to 2 pi
-
-  return angle/ TMath::TwoPi() * nWire;
-
-}
-
 #endif 

TrackRecon.C

@@ -1,6 +1,6 @@
-#define Analyzer_cxx
+#define TrackRecon_cxx
 
-#include "Analyzer.h"
+#include "TrackRecon.h"
 #include <TH2.h>
 #include <TStyle.h>
 #include <TCanvas.h>
@@ -14,20 +14,6 @@
 
 #include "TVector3.h"
 
-TH2F * hsx3IndexVE;
-TH2F * hqqqIndexVE;
-TH2F * hpcIndexVE;
-
-TH2F * hsx3Coin;
-TH2F * hqqqCoin;
-TH2F * hpcCoin;
-
-TH2F * hqqqPolar;
-TH2F * hsx3VpcIndex;
-TH2F * hqqqVpcIndex;
-TH2F * hqqqVpcE;
-TH2F * hsx3VpcE;
-TH2F * hanVScatsum;
 int padID = 0;
 
 SX3 sx3_contr;
@@ -37,7 +23,7 @@ bool HitNonZero;
 
 TH1F * hZProj;
 
-void Analyzer::Begin(TTree * /*tree*/){
+void TrackRecon::Begin(TTree * /*tree*/){
   TString option = GetOption();
 
   hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
@@ -47,17 +33,14 @@ void Analyzer::Begin(TTree * /*tree*/){
 
 }
 
-
+Bool_t TrackRecon::Process(Long64_t entry){
-
-
-Bool_t Analyzer::Process(Long64_t entry){
 
   // if ( entry > 100 ) return kTRUE;
 
   hitPos.Clear();
   HitNonZero = false;
 
-  // if( entry > 1) return kTRUE;
+  if( entry > 1) return kTRUE;
   // printf("################### ev : %llu \n", entry);
 
   b_sx3Multi->GetEntry(entry);
@@ -133,13 +116,6 @@ Bool_t Analyzer::Process(Long64_t entry){
         }else{
           sx3ChBk = sx3.ch[index];
         }
-         for( int j = 0; j < pc.multi; j++){
-      // hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
-            if( sx3.ch[index] > 8 ){
-              hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
-          //  hpcIndexVE->Fill( pc.index[i], pc.e[i] );
-            }
-          }
       }
 
       sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
@@ -170,7 +146,7 @@ Bool_t Analyzer::Process(Long64_t entry){
         double theta = -TMath::Pi()/2 + 2*TMath::Pi()/16/4.*(qqq.id[i]*16 + chWedge +0.5);
         double rho   = 10.+40./16.*(chRing+0.5);
         // if(qqq.e[i]>50){
-        hqqqPolar->Fill( theta, rho);
+        // hqqqPolar->Fill( theta, rho);
         // }
         // qqq.used[i] = true;
         // qqq.used[j] = true;
@@ -187,7 +163,6 @@ Bool_t Analyzer::Process(Long64_t entry){
    
   }
   // //======================= PC
- PCHit_1An hitInfo;
 
   ID.clear();
   int counter=0;
@@ -197,34 +172,52 @@ Bool_t Analyzer::Process(Long64_t entry){
   if( E.size()==3 ){
     float aE = 0;
     float cE = 0;
-    bool multi_an =false;
+    int multi_an =0;
     for(int l=0;l<E.size();l++){
      if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){
-          if(!multi_an){
+        multi_an++;
-            aE = E[l].second;
       }
-          multi_an=true;
-        }
-        else {
-          cE = E[l].second + cE;
-        }
-      }
-    // printf("anode= %d, cathode = %d\n", aID, cID);
-  // }
-    if( ID[0].first < 1 ) {
-      aID = pc.ch[ID[0].second];
-      cID = pc.ch[ID[1].second];
-    }else{
-      cID = pc.ch[ID[0].second];
-      aID = pc.ch[ID[1].second];
     }
 
-    hanVScatsum->Fill(aE,cE);
+    if(multi_an==1){
+      for(int l=0;l<E.size();l++){
+        if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){
+          aE = E[l].second;
+        }else if(E[l].first>24){
+          cE = E[l].second;
+        }
+      }
+    }
+
+  //using CalTrack3 to get the track position and direction
+  
+    // hanVScatsum->Fill(aE,cE);
       
     if( HitNonZero){
-      pw_contr.CalTrack3( hitPos, hitinfo, cID);
+      if (ID.size() == 3) {
+        int aID = -1;
+        int cID1 = -1;
+        int cID2 = -1;
+       for (int i = 0; i < ID.size(); i++) {
+          if (pc.ch[ID[i].second] < 24 && pc.ch[ID[i].second] != 20 && pc.ch[ID[i].second] != 12) {
+          aID = pc.ch[ID[i].second];
+          } else if (pc.ch[ID[i].second] > 24) {
+             if (cID1 == -1) {
+              cID1 = pc.ch[ID[i].second];
+              } else {
+              cID2 = pc.ch[ID[i].second];
+              }
+            }
+         }
+        if (aID != -1 && cID1 != -1 && cID2 != -1) {
+          pw_contr.CalTrack3(hitPos, aID, cID1, cID2);
+          pw_contr.Print();
+          printf("###################\n");
+
           hZProj->Fill(pw_contr.GetZ0());
         }
+      }
+    }
 
   // }
   }
@@ -240,64 +233,13 @@ Bool_t Analyzer::Process(Long64_t entry){
   return kTRUE;
 }
 
-void Analyzer::Terminate(){
+void TrackRecon::Terminate(){
 
   gStyle->SetOptStat("neiou");
-  TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); 
+  TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 200, 200); 
-  canvas->Divide(3,3);
+  padID=1;
-
-  //hsx3VpcIndex->Draw("colz");
-
-  //=============================================== pad-1
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hsx3IndexVE->Draw("colz");
-
-  //=============================================== pad-2
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hqqqIndexVE->Draw("colz");
-
-  //=============================================== pad-3
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hpcIndexVE->Draw("colz");
-
-  //=============================================== pad-4
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hsx3Coin->Draw("colz");
-
-  //=============================================== pad-5
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hqqqCoin->Draw("colz");
-
-  //=============================================== pad-6
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hpcCoin->Draw("colz");
-
-  //=============================================== pad-7
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hsx3VpcIndex ->Draw("colz"); 
-  // hsx3VpcE->Draw("colz") ;
-
-  //=============================================== pad-8
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hqqqVpcIndex ->Draw("colz");  
-
-  // hqqqVpcE ->Draw("colz");
-  //=============================================== pad-9
-  padID ++; 
-
-  // canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
-  // hqqqPolar->Draw("same colz pol");
-
  canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
+
  hZProj->Draw();
-  // hanVScatsum->Draw("colz");
 
 }

TrackRecon.h

@@ -0,0 +1,114 @@
+#ifndef TrackRecon_h
+#define TrackRecon_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TSelector.h>
+
+#include "Armory/ClassDet.h"
+
+class TrackRecon : public TSelector {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   Det sx3;
+   Det qqq;
+   Det pc ;
+
+   ULong64_t       evID;
+   UInt_t          run;
+
+   // List of branches
+   TBranch        *b_eventID;   //!
+   TBranch        *b_run;   //!
+   TBranch        *b_sx3Multi;   //!
+   TBranch        *b_sx3ID;   //!
+   TBranch        *b_sx3Ch;   //!
+   TBranch        *b_sx3E;   //!
+   TBranch        *b_sx3T;   //!
+   TBranch        *b_qqqMulti;   //!
+   TBranch        *b_qqqID;   //!
+   TBranch        *b_qqqCh;   //!
+   TBranch        *b_qqqE;   //!
+   TBranch        *b_qqqT;   //!
+   TBranch        *b_pcMulti;   //!
+   TBranch        *b_pcID;   //!
+   TBranch        *b_pcCh;   //!
+   TBranch        *b_pcE;   //!
+   TBranch        *b_pcT;   //!
+
+   TrackRecon(TTree * /*tree*/ =0) : fChain(0) { }
+   virtual ~TrackRecon() { }
+   virtual Int_t   Version() const { return 2; }
+   virtual void    Begin(TTree *tree);
+   virtual void    SlaveBegin(TTree *tree);
+   virtual void    Init(TTree *tree);
+   virtual Bool_t  Notify();
+   virtual Bool_t  Process(Long64_t entry);
+   virtual Int_t   GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
+   virtual void    SetOption(const char *option) { fOption = option; }
+   virtual void    SetObject(TObject *obj) { fObject = obj; }
+   virtual void    SetInputList(TList *input) { fInput = input; }
+   virtual TList  *GetOutputList() const { return fOutput; }
+   virtual void    SlaveTerminate();
+   virtual void    Terminate();
+
+   ClassDef(TrackRecon,0);
+};
+
+#endif
+
+#ifdef TrackRecon_cxx
+void TrackRecon::Init(TTree *tree){
+
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("evID", &evID, &b_eventID);
+   fChain->SetBranchAddress("run", &run, &b_run);
+
+   sx3.SetDetDimension(24,12);
+   qqq.SetDetDimension(4,32);
+   pc.SetDetDimension(2,24);
+
+   fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
+   fChain->SetBranchAddress("sx3ID",    &sx3.id, &b_sx3ID);
+   fChain->SetBranchAddress("sx3Ch",    &sx3.ch, &b_sx3Ch);
+   fChain->SetBranchAddress("sx3E",     &sx3.e, &b_sx3E);
+   fChain->SetBranchAddress("sx3T",     &sx3.t, &b_sx3T);
+   fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
+   fChain->SetBranchAddress("qqqID",    &qqq.id, &b_qqqID);
+   fChain->SetBranchAddress("qqqCh",    &qqq.ch, &b_qqqCh);
+   fChain->SetBranchAddress("qqqE",     &qqq.e, &b_qqqE);
+   fChain->SetBranchAddress("qqqT",     &qqq.t, &b_qqqT);
+   fChain->SetBranchAddress("pcMulti",  &pc.multi, &b_pcMulti);
+   fChain->SetBranchAddress("pcID",     &pc.id, &b_pcID);
+   fChain->SetBranchAddress("pcCh",     &pc.ch, &b_pcCh);
+   fChain->SetBranchAddress("pcE",      &pc.e, &b_pcE);
+   fChain->SetBranchAddress("pcT",      &pc.t, &b_pcT);
+
+}
+
+Bool_t TrackRecon::Notify(){
+
+   return kTRUE;
+}
+
+void TrackRecon::SlaveBegin(TTree * /*tree*/){
+
+   TString option = GetOption();
+
+}
+
+void TrackRecon::SlaveTerminate(){
+
+}
+
+
+#endif // #ifdef TrackRecon_cxx

gainmatch.C

@@ -0,0 +1,541 @@
+#define gainmatch_cxx
+
+#include "gainmatch.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TCutG.h>
+
+
+#include <utility>
+#include <algorithm>
+
+#include "Armory/ClassSX3.h"
+#include "Armory/ClassPW.h"
+
+#include "TVector3.h"
+
+TH2F * hsx3IndexVE;
+TH2F * hqqqIndexVE;
+TH2F * hqqqIndexVE_cut;
+TH2F * hpcIndexVE;
+
+TH2F * hsx3Coin;
+TH2F * hqqqCoin;
+TH2F * hpcCoin;
+TH2F * hpcCoin_cut;
+
+TH2F * hGoodQQQ;
+TH2F * hGoodQQQRingVWedge;
+
+TH2F * hqqqPolar;
+TH2F * hsx3VpcIndex;
+TH2F * hqqqVpcIndex;
+TH2F * hqqqVpcIndex_cut;
+TH2F * hqqqVpcE;
+TH2F * hqqqVpcE_cut;
+TH2F * hqqqVpcE_cut1;
+TH2F * hqqqVpcE_cut2;
+TH2F * hqqqVpcE_cutCoinc;
+TH2F * hsx3VpcE;
+TH2F * hanVScatsum;
+TH2F * hanVScatsum_cut;
+TH2F * hanVScatsum_cut1;
+TH2F * hanVScatsum_cut2;
+TH2F * hsx3Vsx3;
+TH2F * hsx3uVsx3d_01;
+TH2F * hsx3uVsx3d_23;
+TH2F * hsx3uVsx3d_45;
+TH2F * hsx3uVsx3d_67;
+TH2F * hVCID;
+TH1F *hsx3bk_9_shifted ;
+TH1F *hsx3bk_10_shifted ;
+TH1F *hsx3bk_11_shifted ;
+
+int padID = 0;
+
+TCutG *Coinc_cut_set1;
+//TCutG *crap_cut;
+TCutG *AnCathCoinc_cut;
+TCutG *AnCathCoinc_cut1;
+TCutG *AnCathCoinc_cut2;
+
+SX3 sx3_contr;
+PW pw_contr;
+TVector3 hitPos;
+bool HitNonZero;
+bool inCut;
+bool inCut1;
+bool inCut2;
+bool inCutCoinc;
+TH1F *hZd_01_1;
+TH1F *hZd_01_2;
+TH1F *hZd_01_3;
+TH1F *hZd_01_4;
+TH1F * hZProj;
+TH1F * hsx3bk_11;
+TH1F * hsx3bk_10;
+TH1F * hsx3bk_9;
+TH1F * hsx3bk_8;
+void gainmatch::Begin(TTree * /*tree*/){
+  TString option = GetOption();
+
+  hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy",  24*12, 0,  24*12, 400, 0, 5000); hsx3IndexVE->SetNdivisions( -612, "x");
+  hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
+    hqqqIndexVE_cut = new TH2F("hqqqIndexVE_cut", "QQQ index vs Energy gated; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
+  hpcIndexVE = new TH2F("hpcIndexVE",   "PC index vs Energy; PC index ; Energy",     2*24, 0,   2*24, 400, 0, 4000); hpcIndexVE->SetNdivisions( -1204, "x");
+
+  hGoodQQQ = new TH2F("hGoodQQQ", "number of good QQQ vs QQQ id", 10, 0, 10, 4, 0, 4);
+  hGoodQQQRingVWedge = new TH2F("hGoodQQQRingVWedge", "Ring index, Wedge index", 16*4, 0, 16*4, 16*4, 0, 16*4);
+  hZd_01_1 =new TH1F("hZd_01_1", "Z position", 100, -1, 1); 
+  hZd_01_2 =new TH1F("hZd_01_2", "Z position", 100, -1, 1); 
+  hZd_01_3 =new TH1F("hZd_01_3", "Z position", 100, -1, 1); 
+  hZd_01_4 =new TH1F("hZd_01_4", "Z position", 100, -1, 1); 
+  hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident",  24*12, 0,  24*12,  24*12, 0, 24*12);
+  hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4*2*16, 0, 4*2*16, 4*2*16, 0, 4*2*16);
+  hpcCoin  = new TH2F("hpcCoin",  "PC Coincident",    2*24, 0,   2*24,   2*24, 0, 2*24);
+  hpcCoin_cut  = new TH2F("hpcCoin_cut",  "PC Coincident gated",    2*24, 0,   2*24,   2*24, 0, 2*24);
+    
+  hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16*4, -TMath::Pi(), TMath::Pi(),16, 10, 50);
+
+  hsx3VpcIndex = new TH2F("hsx3Vpcindex", "sx3 vs pc; sx3 index; pc index",  24*12, 0, 24*12, 48, 0, 48);
+  hsx3Vsx3 = new TH2F("hsx3Vsx3", "sx3 vs sx3; sx3 E; sx3 E",  8000, 0, 16000, 8000, 0, 16000);
+  hsx3uVsx3d_01 = new TH2F("hsx3uVsx3d_01", "sx3u vs sx3d; sx3u E; sx3d E",  100, 0, 1, 100, 0, 1);
+  hsx3uVsx3d_23 = new TH2F("hsx3uVsx3d_23", "sx3u vs sx3d; sx3u E; sx3d E",  100, 0, 1, 100, 0, 1);
+  hsx3uVsx3d_45 = new TH2F("hsx3uVsx3d_45", "sx3u vs sx3d; sx3u E; sx3d E",  1000, 0, 5000, 1000, 0, 5000);
+  hsx3uVsx3d_67 = new TH2F("hsx3uVsx3d_67", "sx3u vs sx3d; sx3u E; sx3d E",  1000, 0, 5000, 1000, 0, 5000);
+  hsx3VpcIndex->SetNdivisions( -612, "x");
+  hsx3VpcIndex->SetNdivisions( -12, "y");
+  
+  hqqqVpcIndex = new TH2F("hqqqVpcindex", "qqq vs pc; qqq index; pc index",  4*2*16, 0, 4*2*16, 48, 0, 48);
+  hqqqVpcIndex->SetNdivisions( -612, "x");
+  hqqqVpcIndex->SetNdivisions( -12, "y");
+
+  hqqqVpcIndex_cut = new TH2F("hqqqVpcindex_cut", "qqq vs pc gated; qqq index; pc index",  4*2*16, 0, 4*2*16, 48, 0, 48);
+  hqqqVpcIndex_cut->SetNdivisions( -612, "x");
+  hqqqVpcIndex_cut->SetNdivisions( -12, "y");
+
+  hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
+  hqqqVpcE->SetNdivisions( -612, "x");
+  hqqqVpcE->SetNdivisions( -12, "y");
+    
+  hqqqVpcE_cut = new TH2F("hqqqVpcEnergy_cut", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
+  hqqqVpcE_cut->SetNdivisions( -612, "x");
+  hqqqVpcE_cut->SetNdivisions( -12, "y");
+    
+    hqqqVpcE_cut1 = new TH2F("hqqqVpcEnergy_cut1", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
+    hqqqVpcE_cut1->SetNdivisions( -612, "x");
+    hqqqVpcE_cut1->SetNdivisions( -12, "y");
+      
+    hqqqVpcE_cut2 = new TH2F("hqqqVpcEnergy_cut2", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
+    hqqqVpcE_cut2->SetNdivisions( -612, "x");
+    hqqqVpcE_cut2->SetNdivisions( -12, "y");
+    
+    hqqqVpcE_cutCoinc = new TH2F("hqqqVpcEnergy_cutCoinc", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
+    hqqqVpcE_cutCoinc->SetNdivisions( -612, "x");
+    hqqqVpcE_cutCoinc->SetNdivisions( -12, "y");
+  hsx3bk_8=new TH1F("hsx3bk_8", "hsx3bk_8",1000, 0,5000);
+  hsx3bk_9=new TH1F("hsx3bk_9", "hsx3bk_9",1000, 0,5000);
+  hsx3bk_10=new TH1F("hsx3bk_10", "hsx3bk_10",1000, 0,5000); 
+  hsx3bk_11=new TH1F("hsx3bk_11", "hsx3bk_11",1000, 0,5000);
+  hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000);
+  hsx3VpcE->SetNdivisions( -612, "x");
+  hsx3VpcE->SetNdivisions( -12, "y");
+  hsx3bk_9_shifted  = new TH1F("hsx3bk_9_shifted", "hsx3bk_9",1000, 0,5000);
+  hsx3bk_10_shifted  = new TH1F("hsx3bk_10_shifted", "hsx3bk_9",1000, 0,5000);
+  hsx3bk_11_shifted  = new TH1F("hsx3bk_11_shifted", "hsx3bk_9",1000, 0,5000);
+  hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
+
+  hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 8000,0 , 16000, 8000, 0 , 16000);
+  hanVScatsum_cut = new TH2F("hanVScatsum_cut", "Anode vs Cathode Sum gated; Anode E; Cathode E", 1600,0 , 16000, 1600, 0 , 16000);
+    hanVScatsum_cut1 = new TH2F("hanVScatsum_cut1", "Anode vs Cathode Sum gated; Anode E; Cathode E", 1600,0 , 16000, 1600, 0 , 16000);
+    hanVScatsum_cut2 = new TH2F("hanVScatsum_cut2", "Anode vs Cathode Sum gated; Anode E; Cathode E", 1600,0 , 16000, 1600, 0 , 16000);
+
+  hVCID = new TH2F("hVCID", "Virtual Cathod ID vs total Cath. Energy", 200, 0, 24, 200, 0, 10000);
+
+  sx3_contr.ConstructGeo();
+  pw_contr.ConstructGeo();
+
+  TFile *f3 = new TFile("Coinc_cut_set1.root");
+    //TFile *f4 = new TFile("crap_cut.root");
+  TFile *f = new TFile("AnCathCoinc_cut.root");
+    TFile *f1 = new TFile("AnCathCoinc_cut1.root");
+    TFile *f2 = new TFile("AnCathCoinc_cut2.root");
+
+    
+  Coinc_cut_set1 = (TCutG*)f3->Get("Coinc_cut_set1");
+    //crap_cut = (TCutG*)f4->Get("crap_cut");
+  AnCathCoinc_cut = (TCutG*)f->Get("AnCathCoinc_cut");
+    AnCathCoinc_cut1 = (TCutG*)f1->Get("AnCathCoinc_cut1");
+    AnCathCoinc_cut2 = (TCutG*)f2->Get("AnCathCoinc_cut2");
+
+    
+}
+
+
+
+
+Bool_t gainmatch::Process(Long64_t entry){
+
+  // if ( entry > 100 ) return kTRUE;
+
+  hitPos.Clear();
+  HitNonZero = false;
+    inCut = false;
+
+  // if( entry > 1) return kTRUE;
+  // printf("################### ev : %llu \n", entry);
+
+  b_sx3Multi->GetEntry(entry);
+  b_sx3ID->GetEntry(entry);
+  b_sx3Ch->GetEntry(entry);
+  b_sx3E->GetEntry(entry);
+  b_sx3T->GetEntry(entry);
+  b_qqqMulti->GetEntry(entry);
+  b_qqqID->GetEntry(entry);
+  b_qqqCh->GetEntry(entry);
+  b_qqqE->GetEntry(entry);
+  b_qqqT->GetEntry(entry);
+  b_pcMulti->GetEntry(entry);
+  b_pcID->GetEntry(entry);
+  b_pcCh->GetEntry(entry);
+  b_pcE->GetEntry(entry);
+  b_pcT->GetEntry(entry);
+
+  sx3.CalIndex();
+  qqq.CalIndex();
+  pc.CalIndex();
+
+  // sx3.Print();
+
+  //########################################################### Raw data
+  // //======================= SX3
+
+std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
+for (int i = 0; i < sx3.multi; i++) {
+    ID.push_back(std::pair<int, int>(sx3.id[i], i));
+
+    hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
+
+    for (int j = i + 1; j < sx3.multi; j++) {
+        hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
+    }
+
+    for (int j = 0; j < pc.multi; j++) {
+        hsx3VpcIndex->Fill(sx3.index[i], pc.index[j]);
+    }
+}
+
+if (ID.size() > 0) {
+    std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b) {
+        return a.first < b.first;
+    });
+
+    std::vector<std::pair<int, int>> sx3ID;
+    sx3ID.push_back(ID[0]);
+    bool found = false;
+    for (size_t i = 1; i < ID.size(); i++) {
+        if (ID[i].first == sx3ID.back().first) {
+            sx3ID.push_back(ID[i]);
+            if (sx3ID.size() >= 3) {
+                found = true;
+            }
+        } else {
+            if (!found) {
+                sx3ID.clear();
+                sx3ID.push_back(ID[i]);
+            }
+        }
+    }
+
+    if (found) {
+        int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1;
+        float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0;
+
+        for (size_t i = 0; i < sx3ID.size(); i++) {
+            int index = sx3ID[i].second;
+
+            if (sx3.ch[index] < 8) {
+                if (sx3.ch[index] % 2 == 0) {
+                    sx3ChDn = sx3.ch[index];
+                    sx3EDn = sx3.e[index];
+                } else {
+                    sx3ChUp = sx3.ch[index];
+                    sx3EUp = sx3.e[index];
+                }
+            } else {
+                sx3ChBk = sx3.ch[index];
+                sx3EBk = sx3.e[index];
+            }
+
+            int ch = sx3.ch[index];
+            float energy = sx3.e[index];
+            if (sx3ID[0].first == 9) {
+               float peak8 = 0.0;
+    		float peak9 = 0.0;
+    		int peak10 = 0.0;
+   		float peak11 = 0.0;
+   		float shift9 =0.0;
+   		float shift10 =0.0;
+   		float shift11 =0.0;
+   		int minBin_8 = hsx3bk_8->FindBin(1);
+                int maxBin_8 = hsx3bk_8->FindBin(5000);
+                int maxRangeBinContent_8 = -1;
+                double maxBinCenter_8 = 0.0;
+                int minBin_9 = hsx3bk_9->FindBin(1);
+                int maxBin_9 = hsx3bk_9->FindBin(5000);
+                int maxRangeBinContent_9 = -1;
+                double maxBinCenter_9 = 0.0;
+                int minBin_10 = hsx3bk_10->FindBin(1);
+                int maxBin_10 = hsx3bk_10->FindBin(5000);
+                int maxRangeBinContent_10 = -1;
+                double maxBinCenter_10 = 0.0;
+                int minBin_11 = hsx3bk_11->FindBin(1);
+                int maxBin_11 = hsx3bk_11->FindBin(5000);
+                int maxRangeBinContent_11 = -1;
+                double maxBinCenter_11 = 0.0;
+                if (sx3ChBk == 8) {
+                   
+                   hsx3bk_8->Fill(sx3EBk);
+                    
+
+                    for (int bin = minBin_8; bin <= maxBin_8; ++bin) {
+                       if (hsx3bk_8->GetBinContent(bin) > maxRangeBinContent_8) {
+                         maxRangeBinContent_8 = hsx3bk_8->GetBinContent(bin);
+                         maxBinCenter_8 = hsx3bk_8->GetBinCenter(bin);
+                       }
+                    }
+
+                   //peak8 = hsx3bk_8->GetMaximumBin();
+                   //peak8 = hsx3bk_8->GetMaximumBin();
+                   //printf("peak8: %f\n",  maxBinCenter_8);
+                }
+               //printf("peak8_mm: %f\n",  maxBinCenter);
+               else if (sx3ChBk == 9) {
+                  
+                   hsx3bk_9->Fill(sx3EBk);
+                   for (int bin = minBin_9; bin <= maxBin_9; ++bin) {
+                       if (hsx3bk_9->GetBinContent(bin) > maxRangeBinContent_9) {
+                         maxRangeBinContent_9 = hsx3bk_9->GetBinContent(bin);
+                         maxBinCenter_9 = hsx3bk_9->GetBinCenter(bin);
+                       }
+                    }
+
+                   //peak8 = hsx3bk_8->GetMaximumBin();
+                   //peak8 = hsx3bk_8->GetMaximumBin();
+                    //printf("peak9: %f\n",  maxBinCenter_9);
+                   //hsx3bk_9_shifted->Fill(sx3EBk*0.76);
+                    peak9 = 2097.5/maxBinCenter_9;
+                    //printf("peak9_shift: %f\n",  peak9);
+                    hsx3bk_9_shifted->Fill(sx3EBk*(2097.5/maxBinCenter_9));
+                   //printf("peak9 %d\n", peak9);
+                }
+               else if(sx3ChBk == 10) {
+                    
+                   hsx3bk_10->Fill(sx3EBk);
+                   for (int bin = minBin_10; bin <= maxBin_10; ++bin) {
+                       if (hsx3bk_10->GetBinContent(bin) > maxRangeBinContent_10) {
+                         maxRangeBinContent_10 = hsx3bk_10->GetBinContent(bin);
+                         maxBinCenter_10 = hsx3bk_10->GetBinCenter(bin);
+                       }
+                    }
+
+                   //peak8 = hsx3bk_8->GetMaximumBin();
+                   //peak8 = hsx3bk_8->GetMaximumBin();
+                    //printf("peak10: %f\n",  maxBinCenter_10);
+                   //hsx3bk_9_shifted->Fill(sx3EBk*0.76);
+                    peak10= 2097.5/maxBinCenter_10;
+                    //printf("peak10_shift: %f\n", 1787.5/maxBinCenter_10);
+                    hsx3bk_10_shifted->Fill(sx3EBk*(2097.5/maxBinCenter_10));
+                   //printf("peak9 %d\n", peak9);
+                }
+                   //peak10 = hsx3bk_10->GetMaximumBin();
+                   // printf("peak10 %d\n" ,peak10);
+               
+              else if(sx3ChBk == 11) {
+                   
+                   hsx3bk_11->Fill(sx3EBk);
+                   for (int bin = minBin_11; bin <= maxBin_11; ++bin) {
+                       if (hsx3bk_11->GetBinContent(bin) > maxRangeBinContent_11) {
+                         maxRangeBinContent_11 = hsx3bk_11->GetBinContent(bin);
+                         maxBinCenter_11 = hsx3bk_11->GetBinCenter(bin);
+                       }
+                    }
+
+                   //peak8 = hsx3bk_8->GetMaximumBin();
+                   //peak8 = hsx3bk_8->GetMaximumBin();
+                    //printf("peak9: %f\n",  maxBinCenter_11);
+                   //hsx3bk_9_shifted->Fill(sx3EBk*0.76);
+                    peak11 = 2097.5/maxBinCenter_11;
+                    //printf("peak11_shift: %f\n",  peak11);
+                    hsx3bk_11_shifted->Fill(sx3EBk*(2097.5/maxBinCenter_11));
+                   //printf("peak9 %d\n", peak9);
+                
+                }
+                 
+    	          
+                  float sx3EBk_shifted = 0.0;
+                  float sx3E_u_matched_01 = 0.0;
+                  float sx3E_d_matched_01 = 0.0;
+                  float sx3E_fb_matched_01 = 0.0;
+                  float sx3E_fbu_matched_01 = 0.0;
+                  float sx3E_fbd_matched_01 = 0.0;  
+                  float diff =0.0;
+                  float ratio = 0.0;
+                  float coeff = 0.0;
+                     if (sx3ChBk == 9) {
+                        sx3EBk_shifted = (sx3EBk *(2097.5/maxBinCenter_9));
+                     } else if (sx3ChBk == 10) {
+                        sx3EBk_shifted = (sx3EBk * (2097.5/maxBinCenter_10));
+                     } else if (sx3ChBk == 11) {
+                       sx3EBk_shifted = (sx3EBk * (2097.5/maxBinCenter_11)) ;
+                     } else {
+                     sx3EBk_shifted = sx3EBk; // Use unshifted value for sx3ChBk == 8
+                     }
+                     if ((sx3ChUp == 1 && sx3ChDn == 0)) {
+                        sx3E_u_matched_01= (sx3EUp-0.898729)/0.836243;
+                        //sx3E_u_matched_01= (0.836243*sx3EDn)+0.898729;
+                        sx3E_d_matched_01= (sx3EDn-0.898729)/0.836243;
+                        sx3E_fb_matched_01=(sx3EBk_shifted+9.2423)/0.924773 ;
+                        sx3E_fbu_matched_01=(sx3E_u_matched_01+9.2423)/0.924773 ;
+                        sx3E_fbd_matched_01=(sx3E_d_matched_01+9.2423)/0.924773 ;
+                        diff = sx3E_fb_matched_01 - (sx3EUp+sx3E_fbd_matched_01);
+                        ratio = sx3EUp/sx3E_fbd_matched_01;
+                        coeff = ((sx3EUp+diff) - (sx3E_fbd_matched_01*ratio))/(diff*(1+ratio));
+                     }
+                     
+                 //TH2F *hsx3uVsx3d_01 = nullptr;
+               if (sx3ChBk >=8) { 
+                    
+                  //if (sx3ChBk == 9) {
+                     
+
+                       
+                          if ((sx3ChUp == 1 && sx3ChDn == 0)) {
+                             if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
+                                if (sx3EBk_shifted > 50 && sx3EUp > 50 && sx3EDn > 50) {
+                                 printf("sx3EUp: %f, sx3EDn: %f, sx3E_u_matched_01: %f,sx3E_d_matched_01: %f\n", sx3EUp, sx3EDn, sx3E_u_matched_01,sx3E_d_matched_01);
+                              //printf("Filling hsx3uVsx3d_01_shifted: %f\n", sx3EBk_ud_matched_01 / sx3EBk_shifted);   
+                             // hsx3uVsx3d_01->Fill(sx3E_u_matched_01 / sx3EBk_shifted, sx3E_d_matched_01 / sx3EBk_shifted);
+                                  hsx3uVsx3d_01->Fill(sx3EUp / sx3EBk_shifted, sx3E_d_matched_01 / sx3EBk_shifted);                              
+                                  hsx3uVsx3d_23->Fill(sx3EUp / sx3EBk_shifted, sx3EDn/ sx3EBk_shifted);
+                            }
+                       } 
+                 //}
+                }     
+                  else if ((sx3ChUp == 3 && sx3ChDn == 2)) {
+                       if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
+                         if (sx3EBk_shifted != 0 && sx3EBk_shifted > 50 && sx3EUp > 50 && sx3EDn > 50) {
+                              printf("sx3EUp: %f, sx3EDn: %f, sx3EBk_shifted: %f\n", sx3EUp, sx3EDn, sx3EBk_shifted);
+                              printf("Filling hsx3uVsx3d_23_shifted: %f\n", sx3EUp / sx3EBk_shifted);
+                             // hsx3uVsx3d_23->Fill(sx3EUp / sx3EBk_shifted, (-0.924754*sx3EUp+0.916671) / sx3EBk_shifted);
+                          
+                         }
+                      }    
+                  }   
+                   /* if ((sx3ChUp == 1 && sx3ChDn == 0)) {
+                        if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
+                          if (sx3EBk != 0 && sx3EBk > 50 && sx3EUp > 50 && sx3EDn > 50) {
+                            printf("sx3EUp: %f, sx3EDn: %f, sx3EBk: %f\n", sx3EUp, sx3EDn, sx3EBk);
+                            printf("Filling hsx3uVsx3d_01: %f\n", sx3EUp / sx3EBk);
+                            hsx3uVsx3d_45->Fill(sx3EUp / sx3EBk, sx3EDn / sx3EBk);
+                          
+                         }
+                      }    
+                  } 
+                  else if ((sx3ChUp == 3 && sx3ChDn == 2)) {
+                       if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
+                         if (sx3EBk != 0 && sx3EBk > 50 && sx3EUp > 50 && sx3EDn > 50) {
+                            printf("sx3EUp: %f, sx3EDn: %f, sx3EBk: %f\n", sx3EUp, sx3EDn, sx3EBk);
+                            printf("Filling hsx3uVsx3d_23: %f\n", sx3EUp / sx3EBk);
+                            hsx3uVsx3d_67->Fill(sx3EUp / sx3EBk, sx3EDn / sx3EBk);
+                           
+                         }
+                      }    
+                }*/
+                 if (sx3ChUp == 1 && sx3ChDn == 0){
+                //if (sx3ChUp == 1 || sx3ChDn == 0 || sx3ChUp == 3 || sx3ChDn == 2 || sx3ChUp == 5 || sx3ChDn == 4 || sx3ChUp == 7 || sx3ChDn == 6) {
+                        if (sx3ChUp != -1  && sx3ChBk !=-1 && sx3ChDn !=-1) {
+                           if (sx3EBk_shifted > 50 && sx3EUp > 50 && sx3EDn>50  &&sx3E_u_matched_01>50 && sx3E_u_matched_01>50) {
+                              //printf("sx3EUp: %f, sx3EDn: %f, sx3E_u_matched_01: %f,sx3E_d_matched_01: %f\n", sx3EUp, sx3EDn, sx3E_u_matched_01,sx3E_d_matched_01);
+                              printf("Filling hsx3uVsx3d_nn: %f, gggggg: %f \n", (sx3EUp+sx3EDn),(sx3E_u_matched_01+sx3E_d_matched_01) );   
+                              hsx3uVsx3d_45->Fill((sx3EUp+sx3E_d_matched_01),sx3EBk_shifted);
+                              hsx3uVsx3d_67->Fill((sx3EUp+sx3E_d_matched_01),sx3E_fb_matched_01);
+                            } 
+                       } 
+                 }
+                  /*if (sx3ChBk > 8) {
+                     if ((sx3ChUp == 7 && sx3ChDn == 6) || 
+                         (sx3ChUp == 5 && sx3ChDn == 4) || 
+                         (sx3ChUp == 3 && sx3ChDn == 2) || 
+                         (sx3ChUp == 1 && sx3ChDn == 0)) {
+                    if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
+                        if (sx3EBk != 0 && sx3EBk > 50 && sx3EUp > 50 && sx3EDn > 50) {
+                            hsx3uVsx3d->Fill(sx3EUp / sx3EBk, sx3EDn / sx3EBk);
+                            hsx3Vsx3->Fill(sx3EUp ,sx3EDn);
+                            printf("sx3EUp: %f | sx3EDn: %f | sx3EBk: %f | sx3ChUp: %d | sx3ChDn: %d | sx3ChBk: %d\n", sx3EUp, sx3EDn, sx3EBk, sx3ChUp, sx3ChDn, sx3ChBk);
+                        }
+                    }
+                  }
+                 }*/
+                 //else {
+                    //printf("sx3EUp\n");
+                //}
+            if (sx3ChUp == 1 && sx3ChDn == 0){
+               
+                if (sx3ChUp != -1  && sx3ChBk !=-1 && sx3ChDn !=-1) {
+                    if (sx3E_d_matched_01> sx3EUp ) {
+                       //printf("hZd_01_1_dn: %f\n", sx3E_d_matched_01);
+                        //printf("hZd_01_1_b: %f\n", sx3E_fb_matched_01);
+                        hZd_01_1->Fill((2*(sx3E_d_matched_01+(coeff*diff))/sx3E_fb_matched_01)-1);
+                        }
+                    else if(sx3EUp> sx3E_d_matched_01) {
+                        //printf("hZd_01_2_sx3EUp: %f\n",sx3EUp );
+                        //printf("hZd_01_2_sx3EDn: %f\n",sx3E_fb_matched_01);
+                        
+                        hZd_01_2->Fill(1-(2*(sx3EUp+(1-coeff)*diff))/sx3E_fb_matched_01);
+                        }  
+                    else if(sx3EUp>0.0  && sx3E_d_matched_01>0.0 && sx3E_d_matched_01>=sx3EUp ) {
+                        hZd_01_3->Fill((2*(sx3E_d_matched_01+ coeff*diff)/sx3E_fb_matched_01)-1);
+                        } 
+                    else if(sx3EUp>0.0  && sx3E_d_matched_01>0.0 && sx3E_d_matched_01<sx3EUp ) {
+                        hZd_01_4->Fill(1-(2*(sx3EUp+ (1-coeff)*diff)/sx3E_fb_matched_01));
+                       } 
+                 }
+            }
+          }
+            for (int j = 0; j < pc.multi; j++) {
+                if (sx3.ch[index] > 8) {
+                    hsx3VpcE->Fill(sx3.e[i], pc.e[j]);
+                }
+            }
+        }
+
+        sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
+        hitPos = sx3_contr.GetHitPos();
+        HitNonZero = true;
+        // hitPos.Print();
+    }
+}
+
+
+
+
+
+  // //======================= PC
+
+
+  //########################################################### Track constrcution
+
+
+  //############################## DO THE KINEMATICS
+
+}
+  return kTRUE;
+}
+
+void gainmatch::Terminate(){
+
+ 
+}

gainmatch.h

@@ -0,0 +1,114 @@
+#ifndef gainmatch_h
+#define gainmatch_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TSelector.h>
+
+#include "Armory/ClassDet.h"
+
+class gainmatch : public TSelector {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   Det sx3;
+   Det qqq;
+   Det pc ;
+
+   ULong64_t       evID;
+   UInt_t          run;
+
+   // List of branches
+   TBranch        *b_eventID;   //!
+   TBranch        *b_run;   //!
+   TBranch        *b_sx3Multi;   //!
+   TBranch        *b_sx3ID;   //!
+   TBranch        *b_sx3Ch;   //!
+   TBranch        *b_sx3E;   //!
+   TBranch        *b_sx3T;   //!
+   TBranch        *b_qqqMulti;   //!
+   TBranch        *b_qqqID;   //!
+   TBranch        *b_qqqCh;   //!
+   TBranch        *b_qqqE;   //!
+   TBranch        *b_qqqT;   //!
+   TBranch        *b_pcMulti;   //!
+   TBranch        *b_pcID;   //!
+   TBranch        *b_pcCh;   //!
+   TBranch        *b_pcE;   //!
+   TBranch        *b_pcT;   //!
+
+   gainmatch(TTree * /*tree*/ =0) : fChain(0) { }
+   virtual ~gainmatch() { }
+   virtual Int_t   Version() const { return 2; }
+   virtual void    Begin(TTree *tree);
+   virtual void    SlaveBegin(TTree *tree);
+   virtual void    Init(TTree *tree);
+   virtual Bool_t  Notify();
+   virtual Bool_t  Process(Long64_t entry);
+   virtual Int_t   GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
+   virtual void    SetOption(const char *option) { fOption = option; }
+   virtual void    SetObject(TObject *obj) { fObject = obj; }
+   virtual void    SetInputList(TList *input) { fInput = input; }
+   virtual TList  *GetOutputList() const { return fOutput; }
+   virtual void    SlaveTerminate();
+   virtual void    Terminate();
+
+   ClassDef(gainmatch,0);
+};
+
+#endif
+
+#ifdef gainmatch_cxx
+void gainmatch::Init(TTree *tree){
+
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("evID", &evID, &b_eventID);
+   fChain->SetBranchAddress("run", &run, &b_run);
+
+   sx3.SetDetDimension(24,12);
+   qqq.SetDetDimension(4,32);
+   pc.SetDetDimension(2,24);
+
+   fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
+   fChain->SetBranchAddress("sx3ID",    &sx3.id, &b_sx3ID);
+   fChain->SetBranchAddress("sx3Ch",    &sx3.ch, &b_sx3Ch);
+   fChain->SetBranchAddress("sx3E",     &sx3.e, &b_sx3E);
+   fChain->SetBranchAddress("sx3T",     &sx3.t, &b_sx3T);
+   fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
+   fChain->SetBranchAddress("qqqID",    &qqq.id, &b_qqqID);
+   fChain->SetBranchAddress("qqqCh",    &qqq.ch, &b_qqqCh);
+   fChain->SetBranchAddress("qqqE",     &qqq.e, &b_qqqE);
+   fChain->SetBranchAddress("qqqT",     &qqq.t, &b_qqqT);
+   fChain->SetBranchAddress("pcMulti",  &pc.multi, &b_pcMulti);
+   fChain->SetBranchAddress("pcID",     &pc.id, &b_pcID);
+   fChain->SetBranchAddress("pcCh",     &pc.ch, &b_pcCh);
+   fChain->SetBranchAddress("pcE",      &pc.e, &b_pcE);
+   fChain->SetBranchAddress("pcT",      &pc.t, &b_pcT);
+
+}
+
+Bool_t gainmatch::Notify(){
+
+   return kTRUE;
+}
+
+void gainmatch::SlaveBegin(TTree * /*tree*/){
+
+   TString option = GetOption();
+
+}
+
+void gainmatch::SlaveTerminate(){
+
+}
+
+
+#endif // #ifdef gainmatch_cxx