diff --git a/Analyzer.C b/Analyzer.C
index 66014c2..893f934 100644
--- a/Analyzer.C
+++ b/Analyzer.C
@@ -72,7 +72,7 @@ void Analyzer::Begin(TTree * /*tree*/){
   hsx3VpcE->SetNdivisions( -612, "x");
   hsx3VpcE->SetNdivisions( -12, "y");
 
-  hZProj = new TH1F("hZProj", "Nos of anodes", 20, 0, 19);
+  hZProj = new TH1F("hZProj", "ZProjection", 600,-600, 600);
   hAnodeHits = new TH2F("hAnodeHits", "Anode vs Anode Energy, Anode ID; Anode E", 24,0 , 23, 400, 0 , 20000);
   hAnodeMultiplicity = new TH1F("hAnodeMultiplicity", "Number of Anodes/Event", 40, 0, 40); 
   hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000);
@@ -270,7 +270,10 @@ Bool_t Analyzer::Process(Long64_t entry){
   int aID = 0;
   int cID = 0;
   int anodeCount = 0;
-
+  float cEMax = 0;
+  int cIDMax = 0;
+  float cEnextMax = 0;
+  int cIDnextMax = 0;
   float aE = 0;
   float cE = 0;
 
@@ -291,8 +294,8 @@ Bool_t Analyzer::Process(Long64_t entry){
       //   hpcCoin->Fill( pc.index[i], pc.index[j]);
       // }
 
-      for (int j=0;j<sx3.multi;j++){
-        if(excludeSX3.find(sx3.index[j]) == excludeSX3.end()){
+      // for (int j=0;j<sx3.multi;j++){
+      //   if(excludeSX3.find(sx3.index[j]) == excludeSX3.end()){
 
         hpcIndexVE->Fill( pc.index[i], pc.e[i] );
         for( int j = i+1; j < pc.multi; j++){
@@ -306,8 +309,8 @@ Bool_t Analyzer::Process(Long64_t entry){
                 cathodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
             }
           // }  
-        }
-      }
+      //   }
+      // }
     // hpcIndexVE->Fill( pc.index[i], pc.e[i] );
       hAnodeMultiplicity->Fill(pc.multi);
 
@@ -325,6 +328,15 @@ Bool_t Analyzer::Process(Long64_t entry){
               for (const auto& cathode : cathodeHits) {
                 int cID = cathode.first;
                 float cE = cathode.second;
+                if(cE>cEMax){
+                  cEMax = cE;
+                  cIDMax = cID;
+                }
+                if(cE>cEnextMax && cE<cEMax){
+                  cEnextMax = cE;
+                  cIDnextMax = cID;
+                }
+
                 cESum += cE;
               }
 
@@ -345,12 +357,12 @@ Bool_t Analyzer::Process(Long64_t entry){
     
 
     if( HitNonZero){
-      pw_contr.CalTrack( hitPos, aID, cID);
+      pw_contr.CalTrack1( hitPos, aID, cIDMax, cIDnextMax, cEMax, cEnextMax,1);
       hZProj->Fill(pw_contr.GetZ0());
     }
 
   
-  
+
   
 
 
diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h
index 2cd6a38..620ee90 100644
--- a/Armory/ClassPW.h
+++ b/Armory/ClassPW.h
@@ -63,6 +63,7 @@ public:
   void ConstructGeo();
   void FindWireID(TVector3 pos, TVector3 direction, bool verbose = false);
   void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
+  void CalTrack1(TVector3 sx3Pos, int anodeID, int cathodeID1, int cathodeID2, float cathodeE1, float cathodeE2, bool verbose);
   void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false);
 
   double CircularMean(std::vector<std::pair<int, double>> wireList);
@@ -235,6 +236,27 @@ inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbo
 
 }
 
+inline void PW::CalTrack1(TVector3 sx3Pos, int anodeID, int cathodeID1, int cathodeID2, float cathodeE1, float cathodeE2, bool verbose){
+
+  trackPos = sx3Pos;
+
+  double q1 = cathodeE1;
+  double q2 = cathodeE2;
+  double fracC = (q1) / (q1 + q2);
+  TVector3 shiftC1 = (Ca[cathodeID2].first - Ca[cathodeID1].first) * fracC;
+  TVector3 shiftC2 = (Ca[cathodeID2].second - Ca[cathodeID1].second) * fracC;
+  TVector3 c1 = Ca[cathodeID1].first + shiftC1;
+  TVector3 c2 = Ca[cathodeID1].second + shiftC2;
+
+  TVector3 n1 = (An[anodeID].first).Cross((sx3Pos).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 void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose){
 
   trackPos = sx3Pos;
diff --git a/PCGainMatch.C b/PCGainMatch.C
index 49d669e..5d4a222 100644
--- a/PCGainMatch.C
+++ b/PCGainMatch.C
@@ -148,8 +148,6 @@ Bool_t PCGainMatch::Process(Long64_t entry){
   std::vector<std::pair<int, double>> cathodeHits={};  
   int aID = 0;
   int cID = 0;
-  int anodeCount = 0;
-  int cathodeCount = 0;
   float aE = 0;
   float cE = 0;
 
@@ -167,10 +165,8 @@ Bool_t PCGainMatch::Process(Long64_t entry){
       float cESum = 0;
       if (pc.index[i] < 24 ) {
             anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
-            anodeCount++;
         } else if (pc.index[i] >= 24) {
             cathodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
-            cathodeCount++;
         }
 
       for(int j=i+1;j<pc.multi;j++){
@@ -182,23 +178,29 @@ Bool_t PCGainMatch::Process(Long64_t entry){
       } 
       if (anodeHits.size()==1 && cathodeHits.size() >= 1) {
           
-          // Accumulate total energy from anode hits
           for (const auto& anode : anodeHits) {
-            if(inPCCut){
+          float cESum = 0;
+          // for(int l=0; l<sx3.multi; l++){
+          //   if (sx3.index[l]==80){
+
               int aID = anode.first;
               float aE = anode.second;
-              aESum += aE; // Sum the anode energy
-              hAnodeHits->Fill(aID, aE);
-            }
-          }
-          
-          // Accumulate total energy from cathode hits
-          for (const auto& cathode : cathodeHits) {
-            if(inPCCut){
-              int cID = cathode.first;
-              float cE = cathode.second;
-              cESum += cE; // Sum the cathode energy
-            }
+              aESum += aE;
+
+              for (const auto& cathode : cathodeHits) {
+                int cID = cathode.first;
+                float cE = cathode.second;
+                // if(cE>cEMax){
+                //   cEMax = cE;
+                //   cIDMax = cID;
+                // }
+                // if(cE>cEnextMax && cE<cEMax){
+                //   cEnextMax = cE;
+                //   cIDnextMax = cID;
+                // }
+
+                cESum += cE;
+              }
           }
           
           inCuth = false;