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",

.vscode/settings.json

@@ -1,7 +1,5 @@
 {
   "files.associations": {
-    "NiceMatStyle.C": "cpp",
-    "*.rmd": "markdown",
     "ryanScript.C": "cpp",
     "ryanSelector.C": "cpp",
     "array": "cpp",
@@ -92,9 +90,6 @@
     "processRun.C": "cpp",
     "TrackRecon.C": "cpp",
     "processRuns.C": "cpp",
-    "Analysis.C": "cpp",
+    "Analysis.C": "cpp"
-    "Analyzer1.C": "cpp",
+  }
-    "PCGainMatch.C": "cpp"
-  },
-  "C_Cpp.clang_format_fallbackStyle": "{BasedonStyle: Google, IndentWidth: 2, ColumnLimit: 0}"
 }

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

@@ -1,16 +1,16 @@
 #define Analyzer_cxx
 
 #include "Analyzer.h"
-#include <TCanvas.h>
-#include <TCutG.h>
 #include <TH2.h>
-#include <TMath.h>
 #include <TStyle.h>
-#include <algorithm>
+#include <TCanvas.h>
-#include <utility>
+#include <TMath.h>
+
+#include <utility>
+#include <algorithm>
 
-#include "Armory/ClassPW.h"
 #include "Armory/ClassSX3.h"
+#include "Armory/ClassPW.h"
 
 #include "TVector3.h"
 
@@ -28,17 +28,6 @@ TH2F *hqqqVpcIndex;
 TH2F * hqqqVpcE;
 TH2F * hsx3VpcE;
 TH2F * hanVScatsum;
-TH2F *hICvsSi;
-TH2F *hAnodeHits;
-TH2F *hSiEvsMCPt;
-TH2F *hRfvsMCPt;
-TH1F *hAnodeHits1d;
-TH1F *hPCMultiplicity;
-TH1F *hRFtime;
-TH1F *hSi;
-TH1F *hSi_gated;
-TH1F *hSiMCPt;
-
 int padID = 0;
 
 SX3 sx3_contr;
@@ -47,25 +36,14 @@ TVector3 hitPos;
 bool HitNonZero;
 
 TH1F * hZProj;
-TCutG *PCCoinc;
-
-TCutG *alpha_cut_up;
-TCutG *alpha_cut_down;
-TCutG *cutg;
-bool inCut;
-bool inCutUp;
-bool inCutDown;
-bool inCutG;
 
 void Analyzer::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 = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy",  24*12, 0,  24*12, 400, 0, 5000); hsx3IndexVE->SetNdivisions( -612, "x");
-  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 = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4 * 2 * 16, 0, 4 * 2 * 16, 400, 0, 5000);
+  hpcIndexVE = new TH2F("hpcIndexVE",   "PC index vs Energy; PC index ; Energy",     2*24, 0,   2*24, 400, 0, 4000); hpcIndexVE->SetNdivisions( -1204, "x");
-  hqqqIndexVE->SetNdivisions(-1204, "x");
+
-  hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2 * 24, 0, 2 * 24, 6400, 0, 30000);
-  hpcIndexVE->SetNdivisions(-1204, "x");
 
   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);
@@ -80,42 +58,25 @@ 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, 6400, 0, 30000);
+  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, 6400, 0, 30000);
+  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", "ZProjection", 600, -600, 600);
+  hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
-  hAnodeHits1d = new TH1F("hAnodeHits1d", "Anode Hits", 24, 0, 24);
-  hAnodeHits = new TH2F("hAnodeHits", "Anode vs Anode Energy, Anode ID; Anode E", 24, 0, 23, 400, 0, 30000);
-  hPCMultiplicity = new TH1F("hPCMultiplicity", "Number of PC/Event", 40, 0, 40);
-  hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 6400, 0, 30000, 6400, 0, 30000);
-  hICvsSi = new TH2F("hICvsSi", "IC vs Si; Si E; IC E", 800, 0, 20000, 400, 0, 8000);
-  hSi = new TH1F("hSi", "Si E", 800, 0, 20000);
-  hSi_gated = new TH1F("hSi_gated", "Si E", 800, 0, 20000);
-  hRFtime = new TH1F("hRFtime", "Rf-MCP time(ns)", 3000, -3000, 3000);
-  hSiEvsMCPt = new TH2F("hSiEsMCPt", "Si E vs MCP time; SI E; MCP time", 800, 0, 20000, 3000, -3000, 3000);
-  hSiMCPt = new TH1F("hSiMCPt", "Si vs MCP time", 1500, -3000, 3000);
-  hRfvsMCPt = new TH2F("hRfvsMCPt", "RF vs MCP time; RF(ns) ; MCP time(ns)", 1000, -2000, 2000, 1000, -2000, 2000);
 
+  hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000);
   sx3_contr.ConstructGeo();
   pw_contr.ConstructGeo();
-  TFile *f1 = new TFile("PCCoinc.root");
-  PCCoinc = (TCutG *)f1->Get("PCCoinc");
-  TFile *f2 = new TFile("alpha_cut_up.root");
-  alpha_cut_up = (TCutG *)f2->Get("alpha_cut_up");
-  TFile *f3 = new TFile("alpha_cut_down.root");
-  alpha_cut_down = (TCutG *)f3->Get("alpha_cut_down");
-  TFile *f4 = new TFile("CUTG.root");
-  cutg = (TCutG *)f4->Get("CUTG");
 
-  // TFile *f1 = new TFile("AnCatSum.root");
-  // AnCatSum= (TCutG*)f1->Get("AnCatSum");
 }
 
+
+
+
 Bool_t Analyzer::Process(Long64_t entry){
 
   // if ( entry > 100 ) return kTRUE;
@@ -141,12 +102,6 @@ Bool_t Analyzer::Process(Long64_t entry) {
   b_pcCh->GetEntry(entry);
   b_pcE->GetEntry(entry);
   b_pcT->GetEntry(entry);
-  b_miscCh->GetEntry(entry);
-  b_miscE->GetEntry(entry);
-  b_miscID->GetEntry(entry);
-  b_miscMulti->GetEntry(entry);
-  b_miscT->GetEntry(entry);
-  b_miscTf->GetEntry(entry);
 
   sx3.CalIndex();
   qqq.CalIndex();
@@ -159,7 +114,6 @@ Bool_t Analyzer::Process(Long64_t entry) {
 
   std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
   for( int i = 0; i < sx3.multi; i ++){
-    if (sx3.e[i] > 50) {
     ID.push_back(std::pair<int, int>(sx3.id[i], i));
 
     hsx3IndexVE->Fill( sx3.index[i], sx3.e[i] );
@@ -168,17 +122,19 @@ Bool_t Analyzer::Process(Long64_t entry) {
       hsx3Coin->Fill( sx3.index[i], sx3.index[j]);
     }
 
-      // for( int j = 0; j < pc.multi; j++){
+    for( int j = 0; j < pc.multi; j++){
-      // hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
+      hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
       //  if( sx3.ch[index] > 8 ){
       //    hsx3VpcE->Fill( sx3.e[i], pc.e[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::sort(ID.begin(), ID.end(),  [](const std::pair<int, int> & a, const std::pair<int, int> & b) {
+      return a.first < b.first;
+    } );
     // printf("##############################\n");
     // for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second );
 
@@ -209,6 +165,7 @@ Bool_t Analyzer::Process(Long64_t entry) {
         int index = sx3ID[i].second;
         // printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]);
 
+
         if( sx3.ch[index] < 8 ){
           if( sx3.ch[index] % 2 == 0) {
             sx3ChDn = sx3.ch[index];
@@ -222,9 +179,8 @@ Bool_t Analyzer::Process(Long64_t entry) {
         }
          for( int j = 0; j < pc.multi; j++){
       // hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
-          if (sx3.ch[i] > 8 && pc.index[j] < 24 && pc.e[j] > 50) {
+            if( sx3.ch[index] > 8 ){
               hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
-            // printf(" sx3 Ch: %d, pc Ch: %d , : %d\n", sx3.index[i], pc.index[j], sx3.t[i] - pc.t[j]);
           //  hpcIndexVE->Fill( pc.index[i], pc.e[i] );
             }
           }
@@ -235,33 +191,34 @@ Bool_t Analyzer::Process(Long64_t entry) {
       HitNonZero = true;
       // hitPos.Print();
     }
+
   }
   
+
   // //======================= QQQ
   for( int i = 0; i < qqq.multi; i ++){
-    if (qqq.e[i] > 50) {
-
     // for( int j = 0; j < pc.multi; j++){
         // if(pc.index[j]==4){
           hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] );
         // }
     // }
     for( int j = 0; j < qqq.multi; j++){
-        if (j == i)
+      if ( j == i ) continue;
-          continue;
       hqqqCoin->Fill( qqq.index[i], qqq.index[j]);
     }
 
-      for (int j = 0; j < pc.multi; j++) {
+
-        if (pc.index[j] < 24 && pc.e[j] > 50) {
+    for( int j = i + 1; j < qqq.multi; j++){
-          hqqqVpcE->Fill(qqq.e[i], pc.e[j]);
+      for( int k = 0; k < pc.multi; k++){
+        if(pc.index[k]<24 && pc.e[k]>50 ){
+          hqqqVpcE->Fill( qqq.e[i], pc.e[k] );
           //  hpcIndexVE->Fill( pc.index[i], pc.e[i] );
                 hqqqVpcIndex->Fill( qqq.index[i], pc.index[j] );
-        }
+
         }
       // }
+      }
       // if( qqq.used[i] == true ) continue;
-      for (int j = i + 1; j < qqq.multi; j++) {
       
       //if( qqq.id[i] == qqq.id[j] && (16 - qqq.ch[i]) * (16 - qqq.ch[j]) < 0  ){ // must be same detector and wedge and ring
       if( qqq.id[i] == qqq.id[j]  ){ // must be same detector 
@@ -294,226 +251,94 @@ Bool_t Analyzer::Process(Long64_t entry) {
         }
       }
     }
-    }
+  
+   
   }
   // //======================= PC
 
-  std::vector<std::pair<int, double>> anodeHits;
+  ID.clear();
-  std::vector<std::pair<int, double>> cathodeHits;
+  int counter=0;
-  int aID = 0;
+  std::vector<std::pair<int, double>> E; 
-  int cID = 0;
+  E.clear();
-  float cEMax = 0;
-  int cIDMax = 0;
-  float cEnextMax = 0;
-  int cIDnextMax = 0;
-  float aE = 0;
-  float cE = 0;
-
-  // Define the excluded SX3 and QQQ channels
-  std::unordered_set<int> excludeSX3 = {34, 35, 36, 37, 61, 62, 67, 73, 74, 75, 76, 77, 78, 79, 80, 93, 97, 100, 103, 108, 109, 110, 111, 112};
-  std::unordered_set<int> excludeQQQ = {0, 17, 109, 110, 111, 112, 113, 119, 127, 128};
-
   for( int i = 0; i < pc.multi; i ++){
-    // for(int j=0; j<pc.multi;j++){
-    //   if(pc.id[j]==0){
-    //     anodeCount++;
-    //     }
-    //   }
 
-    if (pc.e[i] > 100 & pc.multi < 7) {
+    if( pc.e[i] > 100 ) ID.push_back(std::pair<int, int>(pc.id[i], i));
-      // hpcIndexVE->Fill( pc.index[i], pc.e[i] );
+    if( pc.e[i] > 100 ) E.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
-      // for( int j = i+1; j < pc.multi; j++){
-      //   hpcCoin->Fill( pc.index[i], pc.index[j]);
-      // }
 
-      // for (int j=0;j<sx3.multi;j++){
+    hpcIndexVE->Fill( pc.index[i], pc.e[i] );
-      //   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++) {
-        inCut = false;
-        if (PCCoinc->IsInside(pc.index[i], pc.index[j])) {
-          inCut = true;
-        }
-        // hpcCoin->Fill( pc.index[i], pc.index[j]);
-      }
-      // if(pc.e[i]>100){
-      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]));
-      }
-      // }
-      //   }
-      // }
-      // hpcIndexVE->Fill( pc.index[i], pc.e[i] );
-    }
-  }
-  hPCMultiplicity->Fill(pc.multi);
-
-  float aESum = 0;
-  float cESum = 0;
-  if (anodeHits.size() == 1 && cathodeHits.size() >= 1) {
-
-    inCutDown = false;
-    inCutUp = false;
-    for (const auto &anode : anodeHits) {
-
-      // for(int l=0; l<sx3.multi; l++){
-      //   if (sx3.index[l]==80){
-
-      aID = anode.first;
-      aE = anode.second;
-      aESum += aE;
-      // printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
-    }
-
-    for (const auto &cathode : cathodeHits) {
-      cID = cathode.first;
-      cE = cathode.second;
-      cESum += cE;
-
-      if (cE > cEMax) {
-        cEMax = cE;
-        cIDMax = cID;
-      }
-      if (cE > cEnextMax && cE < cEMax) {
-        cEnextMax = cE;
-        cIDnextMax = cID;
-      }
-    }
-
-    if (alpha_cut_down->IsInside(aE, cESum)) {
-      inCutDown = true;
-    }
-    if (alpha_cut_up->IsInside(aE, cESum)) {
-      inCutUp = true;
-    }
-
-    // if (inCutUp)
-    // {
-    for (int i = 0; i < pc.multi; i++) {
     for( int j = i+1; j < pc.multi; j++){
       hpcCoin->Fill( pc.index[i], pc.index[j]);
-        hpcIndexVE->Fill(pc.index[i], pc.e[i]);
+    
-      }
-    }
-    // }
-    // if (inCut) {
-    hanVScatsum->Fill(aE, cESum);
-    hAnodeHits->Fill(aID, aE);
-    hAnodeHits1d->Fill(anodeHits.size());
-    // }
-    // }
     }
 
-  // Miscellaneous channels including the Lollipop IC and Si detectors and hot needle IC
+  }
-  // Misc ch 0,1, 2, 3, 4 in order are the LIC, LSi, HNIC-difference, MCP, and Rf
+  //  for( size_t i = 0; i < E.size(); i++) printf("%zu | %d %d \n", i, E[i].first, E[i].second ); 
-  bool timing = false;
-  inCutG = false;
-  double SiE = 0;
-  double SiT = 0;
-  double MCPt = 0;
-  double MCPE = 0;
-  double Rft = 0;
-  double ICt = 0;
-  double ICe = 0;
-  double SiCFDt = 0;
-  for (int i = 0; i < misc.multi; i++) {
-    // if (misc.ch[i] == 1 && misc.e[i] > 10000 && misc.e[i] < 15000) {
-    // if(misc.e[i] > 7500 && misc.e[i]<15000) hSi->Fill(misc.e[i]);
 
-    if (misc.ch[i] == 1) {
+  if( E.size()>=3 ){
-      // hSi->Fill(misc.e[i]);
+
-      SiE = misc.e[i];
+    int aID = 0;
-      SiT = misc.t[i] + misc.tf[i] * 4. / 1000;
+    int cID = 0;
-      // hSi->Fill(misc.e[i]);
+
-    }
+    float aE = 0;
-    if (misc.ch[i] == 2) {
+    float cE = 0;
-      ICt = misc.t[i] + misc.tf[i] * 4. / 1000;
+    // if( ID[0].first < 1 ) {
-      ICe = misc.e[i];
+    //   aID = pc.ch[ID[0].second];
-      hSi->Fill(misc.e[i]);
+    //   cID = pc.ch[ID[1].second];
-    }
+    // }else{
-    if (misc.ch[i] == 3) {
+    //   cID = pc.ch[ID[0].second];
-      // only analyze the first MCP in any event
+    //   aID = pc.ch[ID[1].second];
-      if (MCPt == 0) {
+    // }
-        MCPt = misc.t[i] + misc.tf[i] * 4. / 1000;
+    // printf("anode= %d, cathode = %d\n", aID, cID);
-        MCPE = misc.e[i];
+
-      }
+   for( int k = 0; k < qqq.multi; k++){
-    }
+    if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){
-    if (misc.ch[i] == 4) {
+    
-      // only analyze the first RF in any event
+      int multi_an =0;
-      if (Rft == 0) {
+      for(int l=0;l<E.size();l++){
-        Rft = misc.t[i] + misc.tf[i] * 4. / 1000;
+      if(E[l].first<24){
-      }
+          multi_an++;
-    }
-    if (misc.ch[i] == 5) {
-      if (SiCFDt == 0) {
-        SiCFDt = misc.t[i] + misc.tf[i] * 4. / 1000;
         }
       }
 
-    // hSiEvsMCPt1->Fill(SiE, Rft-MCPt);
+      if(multi_an>=1){
-    // hSiEvsMCPt->Fill(ICe, MCPt - Rft);
+        for(int l=0;l<E.size();l++){
-    if (MCPt != 0 && Rft != 0) {
+          if(E[l].first<24 && E[l].first!=19 && E[l].first!=12){
-      // if (SiE > 10200 && SiE < 12200) {
+            aE = E[l].second;
-      // hRfvsMCPt->Fill(Rft - ICt, MCPt - ICt);
+          }else if(E[l].first>24){
-
+            cE = E[l].second;
-      hSiMCPt->Fill(MCPt - ICt);
-      // if(misc.ch[i] == 2 && misc.e[i] > 1000 && misc.e[i]<2000)
-      hRFtime->Fill(Rft - ICt);
-      // }
-      // printf("RF time : %lld %lld %lld\n", Rft, MCPt, (MCPt - Rft));
-      // }
           }
-    // inCutG = true;
-    // if (misc.ch[i] == 1) hSi->Fill(misc.e[i]);
-
-    // for (int j = 0; j < qqq.multi; j++) {
-      // if (pc.id[j] == 0) {
-        hRfvsMCPt->Fill(Rft-ICt, MCPt -ICt);
-        hSiEvsMCPt->Fill(ICe, MCPt - ICt);
-      // }
-    // }
-    for (int j = i + 1; j < misc.multi; j++) {
-      //   if (cutg->IsInside(misc.e[i], misc.e[j])) {
-      //     inCutG = true;
-      //   })
-
-      if (misc.ch[j] == 4 && misc.ch[i] == 3) {
-        // hRFtime->Fill(misc.t[j]*1. + misc.tf[j] * 4. / 1000 - (misc.t[i]*1. + misc.tf[i] * 4. / 1000));
-
-        if (misc.t[j] + misc.tf[j] * 4. / 1000 - (misc.t[i] + misc.tf[i] * 4. / 1000) > 20 && misc.t[j] + misc.tf[j] * 4. / 1000 - (misc.t[i] + misc.tf[i] * 4. / 1000) < 100) {
-          timing = true;
-        }
-        // printf("RF time : %lld %lld %lld %lld %lld\n", misc.t[i], misc.t[j], misc.tf[i], misc.tf[j], (misc.t[j]*1000 + misc.tf[j]*4 - (misc.t[i]*1000 + misc.tf[i]*4)));
         }
       }
     
-    // for (int j = i + 1; j < misc.multi; j++) {
-    if (timing == true) {
-      // hICvsSi->Fill(misc.e[i], misc.e[j]);
-      if (misc.ch[i] == 1) {
-        hSi_gated->Fill(misc.e[i]);
-        // }
      }
-      // }
    }
+    hanVScatsum->Fill(aE,cE);
+      
+    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];
     }
 
     if( HitNonZero){
-    // pw_contr.CalTrack1( hitPos, aID, cIDMax, cIDnextMax, cEMax, cEnextMax,1);
+      pw_contr.CalTrack( hitPos, aID, cID);
-
-    // pw_contr.CalTrack(hitPos, aID, cID);
       hZProj->Fill(pw_contr.GetZ0());
     }
 
+  
+  }
+  
+
+
   //########################################################### Track constrcution
 
+
   //############################## DO THE KINEMATICS
 
+
   return kTRUE;
 }
 
@@ -521,69 +346,50 @@ void Analyzer::Terminate() {
 
   gStyle->SetOptStat("neiou");
   TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); 
-  // TCanvas *a = new TCanvas("aANASEN", "ANASEN", 800, 600);
   canvas->Divide(3,3);
-  // hRFtime->Draw();
+
-  // TCanvas *b = new TCanvas("bANASEN", "ANASEN", 800, 600);
+  //hsx3VpcIndex->Draw("colz");
-  // // hICvsSi->Draw("colz");
-  // hSi->Draw();
 
   //=============================================== pad-1
-  padID++;
+  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  canvas->cd(padID);
-  canvas->cd(padID)->SetGrid(1);
 
   hsx3IndexVE->Draw("colz");
 
   //=============================================== pad-2
-  padID++;
+  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  canvas->cd(padID);
-  canvas->cd(padID)->SetGrid(1);
 
   hqqqIndexVE->Draw("colz");
 
   //=============================================== pad-3
-  padID++;
+  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  canvas->cd(padID);
-  canvas->cd(padID)->SetGrid(1);
 
   hpcIndexVE->Draw("colz");
 
   //=============================================== pad-4
-  padID++;
+  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  canvas->cd(padID);
-  canvas->cd(padID)->SetGrid(1);
 
   hsx3Coin->Draw("colz");
 
   //=============================================== pad-5
-  padID++;
+  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  canvas->cd(padID);
-  canvas->cd(padID)->SetGrid(1);
 
   canvas->cd(padID)->SetLogz(true);
 
   hqqqCoin->Draw("colz");
 
   //=============================================== pad-6
-  padID++;
+  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  canvas->cd(padID);
-  canvas->cd(padID)->SetGrid(1);
 
   hpcCoin->Draw("colz");
 
   //=============================================== pad-7
-  padID++;
+  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  canvas->cd(padID);
-  canvas->cd(padID)->SetGrid(1);
 
   // hsx3VpcIndex ->Draw("colz"); 
   hsx3VpcE->Draw("colz") ;
 
   //=============================================== pad-8
-  padID++;
+  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  canvas->cd(padID);
-  canvas->cd(padID)->SetGrid(1);
 
   // hqqqVpcIndex ->Draw("colz");  
 
@@ -594,10 +400,8 @@ void Analyzer::Terminate() {
   // canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
   // hqqqPolar->Draw("same colz pol");
 
-  canvas->cd(padID);
+ canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  canvas->cd(padID)->SetGrid(1);
 //  hZProj->Draw();
   hanVScatsum->Draw("colz");
-  // hAnodeHits->Draw("colz");
+
-  // // hAnodeMultiplicity->Draw();
 }

Analyzer.h

@@ -18,7 +18,6 @@ public :
    Det sx3;
    Det qqq;
    Det pc ;
-   Det misc;
 
    ULong64_t       evID;
    UInt_t          run;
@@ -41,13 +40,6 @@ public :
    TBranch        *b_pcCh;   //!
    TBranch        *b_pcE;   //!
    TBranch        *b_pcT;   //!
-   TBranch        *b_miscMulti;   //!
-   TBranch        *b_miscID;   //!
-   TBranch        *b_miscCh;   //!
-   TBranch        *b_miscE;   //!
-   TBranch        *b_miscT;   //!
-   TBranch       *b_miscTf;   //!
-
 
    Analyzer(TTree * /*tree*/ =0) : fChain(0) { }
    virtual ~Analyzer() { }
@@ -100,13 +92,6 @@ void Analyzer::Init(TTree *tree){
    fChain->SetBranchAddress("pcCh",     &pc.ch, &b_pcCh);
    fChain->SetBranchAddress("pcE",      &pc.e, &b_pcE);
    fChain->SetBranchAddress("pcT",      &pc.t, &b_pcT);
-   fChain->SetBranchAddress("miscMulti", &misc.multi, &b_miscMulti);
-   fChain->SetBranchAddress("miscID",    &misc.id, &b_miscID);
-   fChain->SetBranchAddress("miscCh",    &misc.ch, &b_miscCh);
-   fChain->SetBranchAddress("miscE",     &misc.e, &b_miscE);
-   fChain->SetBranchAddress("miscT",     &misc.t, &b_miscT);
-   fChain->SetBranchAddress("miscF",     &misc.tf, &b_miscTf);
-
 
 }
 

Analyzer1.C

@@ -1,402 +0,0 @@
-#define Analyzer1_cxx
-
-#include "Analyzer1.h"
-#include <TH2.h>
-#include <TStyle.h>
-#include <TCanvas.h>
-#include <TMath.h>
-
-#include <utility>
-#include <algorithm>
-
-#include "Armory/ClassSX3.h"
-#include "Armory/ClassPW.h"
-
-#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;
-PW pw_contr;
-TVector3 hitPos;
-bool HitNonZero;
-
-TH1F * hZProj;
-
-void Analyzer1::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");
-  hpcIndexVE = new TH2F("hpcIndexVE",   "PC index vs Energy; PC index ; Energy",     2*24, 0,   2*24, 400, 0, 4000); hpcIndexVE->SetNdivisions( -1204, "x");
-
-
-  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);
-
-  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);
-  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");
-
-  hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 400, 0, 5000);
-  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->SetNdivisions( -612, "x");
-  hsx3VpcE->SetNdivisions( -12, "y");
-
-  hZProj = new TH1F("hZProj", "Z Projection", 1200, -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();
-
-}
-
-
-
-
-Bool_t Analyzer1::Process(Long64_t entry){
-
-  // if ( entry > 100 ) return kTRUE;
-
-  hitPos.Clear();
-  HitNonZero = 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( sx3.ch[index] > 8 ){
-      //    hsx3VpcE->Fill( sx3.e[i], pc.e[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;
-    } );
-    // printf("##############################\n");
-    // for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second );
-
-    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]);
-        }
-      }
-    }
-
-    // printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
-
-    if( found ){
-      int sx3ChUp, sx3ChDn, sx3ChBk;
-      float sx3EUp, sx3EDn;
-      // printf("------ sx3 ID : %d, multi: %zu\n", sx3ID[0].first, sx3ID.size());
-      for( size_t i = 0; i < sx3ID.size(); i++ ){
-        int index = sx3ID[i].second;
-        // printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]);
-
-
-        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];
-        }
-         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);
-      hitPos = sx3_contr.GetHitPos();
-      HitNonZero = true;
-      // hitPos.Print();
-    }
-
-  }
-  
-
-  // //======================= QQQ
-  for( int i = 0; i < qqq.multi; i ++){
-    // for( int j = 0; j < pc.multi; j++){
-        // if(pc.index[j]==4){
-          hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] );
-        // }
-    // }
-    for( int j = 0; j < qqq.multi; j++){
-      if ( j == i ) continue;
-      hqqqCoin->Fill( qqq.index[i], qqq.index[j]);
-    }
-
-
-    for( int j = i + 1; j < qqq.multi; j++){
-      for( int k = 0; k < pc.multi; k++){
-        if(pc.index[k]<24 && pc.e[k]>50 ){
-          hqqqVpcE->Fill( qqq.e[i], pc.e[k] );
-          //  hpcIndexVE->Fill( pc.index[i], pc.e[i] );
-                hqqqVpcIndex->Fill( qqq.index[i], pc.index[j] );
-
-        }
-      // }
-      }
-      // if( qqq.used[i] == true ) continue;
-      
-      //if( qqq.id[i] == qqq.id[j] && (16 - qqq.ch[i]) * (16 - qqq.ch[j]) < 0  ){ // must be same detector and wedge and ring
-      if( qqq.id[i] == qqq.id[j]  ){ // must be same detector 
-        
-        int chWedge = -1;
-        int chRing  =  -1;
-        if( qqq.ch[i] < qqq.ch[j]){
-          chRing = qqq.ch[j] - 16;
-          chWedge = qqq.ch[i];
-        }else{
-          chRing = qqq.ch[i];
-          chWedge = qqq.ch[j] - 16;
-        }
-
-        // printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing);
-
-        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);
-        // }
-        // qqq.used[i] = true;
-        // qqq.used[j] = true;
-
-        if( !HitNonZero ){
-        double x = rho * TMath::Cos(theta);
-        double y = rho * TMath::Sin(theta);
-        hitPos.SetXYZ(x, y, 23 + 75 + 30);
-        HitNonZero = true;
-        }
-      }
-    }
-  
-   
-  }
-  // //======================= PC
-
-  ID.clear();
-  int counter=0;
-  std::vector<std::pair<int, double>> E; 
-  E.clear();
-  for( int i = 0; i < pc.multi; i ++){
-
-    if( pc.e[i] > 100 ) ID.push_back(std::pair<int, int>(pc.id[i], i));
-    if( pc.e[i] > 100 ) E.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
-
-    hpcIndexVE->Fill( pc.index[i], pc.e[i] );
-
-    for( int j = i+1; j < pc.multi; j++){
-      hpcCoin->Fill( pc.index[i], pc.index[j]);
-    
-    }
-
-  }
-  //  for( size_t i = 0; i < E.size(); i++) printf("%zu | %d %d \n", i, E[i].first, E[i].second ); 
-
-  if( E.size()>=3 ){
-
-    int aID = 0;
-    int cID = 0;
-
-    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];
-    // }else{
-    //   cID = pc.ch[ID[0].second];
-    //   aID = pc.ch[ID[1].second];
-    // }
-    // printf("anode= %d, cathode = %d\n", aID, cID);
-
-  // for( int k = 0; k < qqq.multi; k++){
-  //   if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){
-      for(int l=0;l<E.size();l++){
-        if(E[l].first<24 ){
-          if(!multi_an){
-            aE = E[l].second;
-          }
-          multi_an=true;
-        }
-        else if (E[l].first>=24){
-          cE = E[l].second + cE;
-        }
-      }
-    // }
-  // }
-    hanVScatsum->Fill(aE,cE);
-      
-    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];
-    }
-
-    if( HitNonZero){
-      pw_contr.CalTrack( hitPos, aID, cID);
-      hZProj->Fill(pw_contr.GetZ0());
-    }
-
-  
-  }
-  
-
-
-  //########################################################### Track constrcution
-
-
-  //############################## DO THE KINEMATICS
-
-
-  return kTRUE;
-}
-
-void Analyzer1::Terminate(){
-
-  gStyle->SetOptStat("neiou");
-  TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); 
-  canvas->Divide(3,3);
-
-  //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);
-
-  canvas->cd(padID)->SetLogz(true);
-
-  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");
-
-}

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/ClassDet.h

@@ -14,7 +14,7 @@ public:
   unsigned short       ch[MAXMULTI];
   unsigned short       e[MAXMULTI];
   unsigned long long   t[MAXMULTI];
-  unsigned long long tf[MAXMULTI];
+
   unsigned short         sn[MAXMULTI];
   unsigned short     digiCh[MAXMULTI];
 
@@ -28,7 +28,6 @@ public:
       ch[i] = 0;
        e[i] = 0;
        t[i] = 0;
-       tf[i] = 0;
        index[i] = 0;
        sn[i] = 0;
        digiCh[i] = 0;

Armory/Mapper.cpp

@@ -67,7 +67,6 @@ int main(int argc, char **argv){
   Det sx3;
   Det qqq;
   Det pc ;
-  Det misc;
 
   printf(" Raw root file : %s\n", inFileName.c_str());
   printf("           Run : %03d\n", run);
@@ -100,14 +99,6 @@ int main(int argc, char **argv){
   newTree->Branch("pcE",     &pc.e,     "pcEnergy[pcMulti]/s");
   newTree->Branch("pcT",     &pc.t,     "pcTime[pcMulti]/l");
 
-  newTree->Branch("miscMulti", &misc.multi, "miscMulti/s");
-  newTree->Branch("miscID",    &misc.id,    "miscID[miscMulti]/s");
-  newTree->Branch("miscCh",    &misc.ch,    "miscCh[miscMulti]/s");
-  newTree->Branch("miscE",     &misc.e,     "miscEnergy[miscMulti]/s");
-  newTree->Branch("miscT",     &misc.t,     "miscTime[miscMulti]/l");
-  newTree->Branch("miscF",     &misc.tf,     "miscFineTime[miscMulti]/l");
-
-
   ///================== looping old tree and apply mapping
 
   //clock   
@@ -121,12 +112,8 @@ int main(int argc, char **argv){
     sx3.multi = 0;
     qqq.multi = 0;
     pc.multi = 0;
-    misc.multi=0;
 
-    sx3.Clear();
     qqq.Clear();
-    pc.Clear();
-    misc.Clear();
 
     for( unsigned int i = 0; i < multi; i++){
 
@@ -173,17 +160,6 @@ int main(int argc, char **argv){
         pc.t[pc.multi] = e_t[i];
         pc.multi ++;
       }
-
-      //=================================== misc
-      if( 30000 <= ID && ID < 40000 ) {
-        misc.id[misc.multi] = (ID - 30000) / 100;
-        misc.ch[misc.multi] = (ID - 30000) % 100;
-        misc.e[misc.multi] = e[i];
-        misc.t[misc.multi] = e_t[i];
-        misc.tf[misc.multi] = e_f[i];
-        // if( ID == 30002 || ID == 30004 ) printf("sn : %d ch: %2d | gID %3d | ID %6d | e_f : %d\n", sn[i], ch[i], globalCh, ID, e_f[i]);
-        misc.multi ++;
-      }
     }
     
     saveFile->cd(); //set focus on this file

PCGainMatch.C

@@ -1,454 +0,0 @@
-#define PCGainMatch_cxx
-
-#include "PCGainMatch.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 * hpcIndexVE;
-
-TH2F * hsx3Coin;
-TH2F * hqqqCoin;
-TH2F * hpcCoin;
-
-TH2F * hqqqPolar;
-TH2F * hsx3VpcIndex;
-TH2F * hqqqVpcIndex;
-TH2F * hqqqVpcE;
-TH2F * hsx3VpcE;
-TH2F * hanVScatsum;
-TH2F * hanVScatsum_a[24];
-TH2F * hanVScatsum_hcut;
-TH2F * hanVScatsum_lcut;
-TH2F * hAnodeHits;
-TH1F * hAnodeMultiplicity;
-
-
-int padID = 0;
-
-SX3 sx3_contr;
-PW pw_contr;
-TVector3 hitPos;
-bool HitNonZero;
-
-TH1F * hZProj;
-TCutG *AnCatSum_high;
-TCutG *AnCatSum_low;
-TCutG *PCCoinc_cut1;
-TCutG *PCCoinc_cut2;
-bool inCuth;
-bool inCutl;
-bool inPCCut;
-
-void PCGainMatch::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");
-  hpcIndexVE = new TH2F("hpcIndexVE",   "PC index vs Energy; PC index ; Energy",     2*24, 0,   2*24, 800, 0, 16000); hpcIndexVE->SetNdivisions( -1204, "x");
-
-
-  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);
-
-  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);
-  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");
-
-  hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 400, 0, 5000);
-  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->SetNdivisions( -612, "x");
-  hsx3VpcE->SetNdivisions( -12, "y");
-
-  hZProj = new TH1F("hZProj", "Nos of anodes", 20, 0, 19);
-  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, 800, 0 , 16000);
-  for (int i = 0; i < 24; i++) {
-        TString histName = Form("hAnodeVsCathode_%d", i);
-        TString histTitle = Form("Anode %d vs Cathode Sum; Anode E; Cathode Sum E", i);
-        hanVScatsum_a[i] = new TH2F(histName, histTitle, 400, 0, 10000, 400, 0, 16000);
-    }
-  hanVScatsum_lcut = new TH2F("hanVScatsum_LCUT", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 16000, 400, 0 , 16000);
-  hanVScatsum_hcut = new TH2F("hanVScatsum_HCUT", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 16000, 400, 0 , 16000);
-
-  sx3_contr.ConstructGeo();
-  pw_contr.ConstructGeo();
-  // TFile *f1 = new TFile("AnCatSum_high.root");
-  // TFile *f2 = new TFile("AnCatSum_low.root");
-  // TFile *f3 = new TFile("PCCoinc_cut1.root");
-  // TFile *f4 = new TFile("PCCoinc_cut2.root");
-  // AnCatSum_high= (TCutG*)f1->Get("AnCatSum_high");
-  // AnCatSum_low= (TCutG*)f2->Get("AnCatSum_low");
-  // PCCoinc_cut1= (TCutG*)f3->Get("PCCoinc_cut1");
-  // PCCoinc_cut2= (TCutG*)f4->Get("PCCoinc_cut2");
-}
-
-
-
-
-Bool_t PCGainMatch::Process(Long64_t entry){
-// if (entry % 1000000 == 0) {
-//     std::cout << "Processing entry: " << entry << std::endl;
-// }
-  // if ( entry > 100 ) return kTRUE;
-
-  hitPos.Clear();
-  HitNonZero = 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
-
-  // //======================= PC
-
-  std::vector<std::pair<int, double>> anodeHits={};
-  std::vector<std::pair<int, double>> cathodeHits={};  
-  int aID = 0;
-  int cID = 0;
-  float aE = 0;
-  float cE = 0;
-
-// Define the excluded SX3 and QQQ channels
-  // std::unordered_set<int> excludeSX3 = {34, 35, 36, 37, 61, 62, 67, 73, 74, 75, 76, 77, 78, 79, 80, 93, 97, 100, 103, 108, 109, 110, 111, 112};
-  // std::unordered_set<int> excludeQQQ = {0, 17, 109, 110, 111, 112, 113, 119, 127, 128};
-  // inCuth=false;
-  // inCutl=false;
-  // inPCCut=false;
-  for( int i = 0; i < pc.multi; i ++){
-    
-    if(pc.e[i]>50 && pc.multi<7){ 
-
-      float aESum = 0;
-      float cESum = 0;
-      if (pc.index[i] < 24 ) {
-            anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
-        } else if (pc.index[i] >= 24) {
-            cathodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
-        }
-
-      for(int j=i+1;j<pc.multi;j++){
-        // if(PCCoinc_cut1->IsInside(pc.index[i], pc.index[j]) || PCCoinc_cut2->IsInside(pc.index[i], pc.index[j])){
-        //   // hpcCoin->Fill(pc.index[i], pc.index[j]);
-        //   inPCCut = true;
-        // }
-        hpcCoin->Fill(pc.index[i], pc.index[j]);
-      } 
-      if (anodeHits.size()==1 && cathodeHits.size() >= 1) {
-          
-          for (const auto& anode : anodeHits) {
-          // for(int l=0; l<sx3.multi; l++){
-          //   if (sx3.index[l]==80){
-
-              aID = anode.first;
-              aE = anode.second;
-              aESum += aE;
-              printf("aID : %d, aE : %f\n", aID, aE);
-          }
-
-            printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
-            for (const auto& cathode : cathodeHits) {
-              cID = cathode.first;
-              cE = cathode.second;
-              // if(cE>cEMax){
-              //   cEMax = cE;
-              //   cIDMax = cID;
-              // }
-              // if(cE>cEnextMax && cE<cEMax){
-              //   cEnextMax = cE;
-              //   cIDnextMax = cID;
-              // }
-
-              cESum += cE;
-              }
-          // }
-          
-          // inCuth = false;
-          // inCutl = false;
-          // inPCCut = false;
-          // for(int j=i+1;j<pc.multi;j++){
-          //   if(PCCoinc_cut1->IsInside(pc.index[i], pc.index[j]) || PCCoinc_cut2->IsInside(pc.index[i], pc.index[j])){
-          //     // hpcCoin->Fill(pc.index[i], pc.index[j]);
-          //     inPCCut = true;
-          //   }
-          //   hpcCoin->Fill(pc.index[i], pc.index[j]);
-          // } 
-
-          // Check if the accumulated energies are within the defined ranges
-          // if (AnCatSum_high && AnCatSum_high->IsInside(aESum, cESum)) {
-          //     inCuth = true;
-          // }
-          // if (AnCatSum_low && AnCatSum_low->IsInside(aESum, cESum)) {
-          //     inCutl = true;
-          // }
-          
-          // Fill histograms based on the cut conditions
-          // if (inCuth && inPCCut) {
-          //     hanVScatsum_hcut->Fill(aESum, cESum);
-          // }
-          // if (inCutl && inPCCut) {
-          //     hanVScatsum_lcut->Fill(aESum, cESum);
-          // }               
-            // for(auto anode : anodeHits){
-         
-              // float aE = anode.second;
-              // aESum += aE;
-              // if(inPCCut){
-            hanVScatsum->Fill(aESum, cESum);
-            // }
-            if (aID < 24 && aE > 50) {
-             hanVScatsum_a[aID]->Fill(aE, cESum);
-            }
-
-          
-            // }
-          // Fill histograms for the `pc` data
-          hpcIndexVE->Fill(pc.index[i], pc.e[i]);
-          // if(inPCCut){
-            hAnodeMultiplicity->Fill(anodeHits.size());
-          // }
-      }
-
-    }
-  }    
-  
-
-
-  // //======================= SX3
-
-  std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
-  for( int i = 0; i < sx3.multi; i ++){
-    if(sx3.e[i]>50){
-      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( sx3.ch[index] > 8 ){
-        //    hsx3VpcE->Fill( sx3.e[i], pc.e[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;
-    } );
-    // printf("##############################\n");
-    // for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second );
-
-    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]);
-        }
-      }
-    }
-
-    // printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
-
-    if( found ){
-      int sx3ChUp, sx3ChDn, sx3ChBk;
-      float sx3EUp, sx3EDn;
-      // printf("------ sx3 ID : %d, multi: %zu\n", sx3ID[0].first, sx3ID.size());
-      for( size_t i = 0; i < sx3ID.size(); i++ ){
-        int index = sx3ID[i].second;
-        // printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]);
-
-
-        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];
-        }
-         for( int j = 0; j < pc.multi; j++){
-
-      // hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
-            if( sx3.ch[index] > 8 && pc.index[j]<24 && pc.e[j]>50 ){
-              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);
-      hitPos = sx3_contr.GetHitPos();
-      HitNonZero = true;
-      // hitPos.Print();
-    }
-
-  }
-  
-
-  // //======================= QQQ
-  for( int i = 0; i < qqq.multi; i ++){
-
-      // for( int j = 0; j < pc.multi; j++){
-        if(qqq.e[i]>50 ){
-          hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] );
-          for( int j = 0; j < qqq.multi; j++){
-            if ( j == i ) continue;
-            hqqqCoin->Fill( qqq.index[i], qqq.index[j]);
-          }
-
-
-          for( int j = i + 1; j < qqq.multi; j++){
-            for( int k = 0; k < pc.multi; k++){
-              // if(qqq.e[i>50]){
-                hqqqVpcE->Fill( qqq.e[i], pc.e[k] );
-                hqqqVpcIndex->Fill( qqq.index[i], pc.index[j] );
-              }
-            // }
-          }
-        }  
-      // }
-  }
-
-    // hanVScatsum->Fill(aE,cE);
-    
-
-    if( HitNonZero){
-      pw_contr.CalTrack( hitPos, aID, cID);
-      hZProj->Fill(pw_contr.GetZ0());
-    }
-
-  
-  
-  
-
-
-  //########################################################### Track constrcution
-
-
-  //############################## DO THE KINEMATICS
-
-
-  return kTRUE;
-}
-
-void PCGainMatch::Terminate(){
-
-  gStyle->SetOptStat("neiou");
-  TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); 
-  canvas->Divide(3,3);
-
-  //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);
-
-  canvas->cd(padID)->SetLogz(true);
-
-  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);
-  hanVScatsum->Draw("colz");
-  // hAnodeHits->Draw("colz");
-  // hAnodeMultiplicity->Draw();
-}

PreAnalysis_C_ACLiC_dict.cxx_tmp_829464

@@ -1,4 +0,0 @@
-// Do NOT change. Changes will be lost next time file is generated
-
-#define R__DICTIONARY_FILENAME dIhomedIvigneshsitaramandI27Al_analysisdIANASEN_analysisdIPreAnalysis_C_ACLiC_dict
-#define R__NO_DEPRECATION

ProcessRun.sh

@@ -1,8 +1,7 @@
 #!/bin/bash
 
-if [ "$#" -ne 3 ]; then
+if [ "$#" -ne 2 ]; then
-    echo "Usage: $0 runID timeWindow_ns option"
+    echo "Usage: $0 runID timeWindow_ns"
-    echo "option: 0 - process raw data, 1 - process mapped data"
     echo "Exiting..."
     exit 1
 fi
@@ -10,24 +9,19 @@ fi
 runID=$1
 timeWindow=$2
 
-option=$3
+rawFolder=/home/tandem/Desktop/analysis/data
+rootFolder=/home/tandem/Desktop/analysis/data/root_data
 
-rawFolder=/home/tandem/data1/2024_09_17Fap/data
+rsync -a splitpole@128.186.111.223:/media/nvmeData/ANASEN27Alap/*.fsu /home/tandem/Desktop/analysis/data
-rootFolder=/home/tandem/data1/2024_09_17Fap/data/root_data
 
-if [ $option -eq 0 ]; then
+fileList=`\ls -1 ${rawFolder}/Run_${runID}_*.fsu`
 
-    rsync -auh --info=progress2 splitpole@128.186.111.223:/media/nvmeData/2024_09_17Fap/*.fsu /home/tandem/data1/2024_09_17Fap/data
+./EventBuilder ${timeWindow} 0 0 10000000 ${fileList}
-
-    fileList=`\ls -1 ${rawFolder}/*Run_${runID}_*.fsu`
-
-    ./EventBuilder ${timeWindow} 0 0 100000000 ${fileList}
 
 outFile=${rawFolder}/*${runID}*${timeWindow}.root
 
 mv -vf ${outFile} ${rootFolder}/.
 
 ./Mapper ${rootFolder}/*${runID}*${timeWindow}.root
-fi
 
 root "processRun.C(\"${rootFolder}/Run_${runID}_mapped.root\")"

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

@@ -1,5 +1,5 @@
-#ifndef Analyzer1_h
+#ifndef TrackRecon_h
-#define Analyzer1_h
+#define TrackRecon_h
 
 #include <TROOT.h>
 #include <TChain.h>
@@ -8,7 +8,7 @@
 
 #include "Armory/ClassDet.h"
 
-class Analyzer1 : public TSelector {
+class TrackRecon : public TSelector {
 public :
    TTree          *fChain;   //!pointer to the analyzed TTree or TChain
 
@@ -41,8 +41,8 @@ public :
    TBranch        *b_pcE;   //!
    TBranch        *b_pcT;   //!
 
-   Analyzer1(TTree * /*tree*/ =0) : fChain(0) { }
+   TrackRecon(TTree * /*tree*/ =0) : fChain(0) { }
-   virtual ~Analyzer1() { }
+   virtual ~TrackRecon() { }
    virtual Int_t   Version() const { return 2; }
    virtual void    Begin(TTree *tree);
    virtual void    SlaveBegin(TTree *tree);
@@ -57,13 +57,13 @@ public :
    virtual void    SlaveTerminate();
    virtual void    Terminate();
 
-   ClassDef(Analyzer1,0);
+   ClassDef(TrackRecon,0);
 };
 
 #endif
 
-#ifdef Analyzer1_cxx
+#ifdef TrackRecon_cxx
-void Analyzer1::Init(TTree *tree){
+void TrackRecon::Init(TTree *tree){
 
    // Set branch addresses and branch pointers
    if (!tree) return;
@@ -95,20 +95,20 @@ void Analyzer1::Init(TTree *tree){
 
 }
 
-Bool_t Analyzer1::Notify(){
+Bool_t TrackRecon::Notify(){
 
    return kTRUE;
 }
 
-void Analyzer1::SlaveBegin(TTree * /*tree*/){
+void TrackRecon::SlaveBegin(TTree * /*tree*/){
 
    TString option = GetOption();
 
 }
 
-void Analyzer1::SlaveTerminate(){
+void TrackRecon::SlaveTerminate(){
 
 }
 
 
-#endif // #ifdef Analyzer_cxx
+#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

@@ -1,5 +1,5 @@
-#ifndef PCGainMatch_h
+#ifndef gainmatch_h
-#define PCGainMatch_h
+#define gainmatch_h
 
 #include <TROOT.h>
 #include <TChain.h>
@@ -8,7 +8,7 @@
 
 #include "Armory/ClassDet.h"
 
-class PCGainMatch : public TSelector {
+class gainmatch : public TSelector {
 public :
    TTree          *fChain;   //!pointer to the analyzed TTree or TChain
 
@@ -41,8 +41,8 @@ public :
    TBranch        *b_pcE;   //!
    TBranch        *b_pcT;   //!
 
-   PCGainMatch(TTree * /*tree*/ =0) : fChain(0) { }
+   gainmatch(TTree * /*tree*/ =0) : fChain(0) { }
-   virtual ~PCGainMatch() { }
+   virtual ~gainmatch() { }
    virtual Int_t   Version() const { return 2; }
    virtual void    Begin(TTree *tree);
    virtual void    SlaveBegin(TTree *tree);
@@ -57,13 +57,13 @@ public :
    virtual void    SlaveTerminate();
    virtual void    Terminate();
 
-   ClassDef(PCGainMatch,0);
+   ClassDef(gainmatch,0);
 };
 
 #endif
 
-#ifdef PCGainMatch_cxx
+#ifdef gainmatch_cxx
-void PCGainMatch::Init(TTree *tree){
+void gainmatch::Init(TTree *tree){
 
    // Set branch addresses and branch pointers
    if (!tree) return;
@@ -95,20 +95,20 @@ void PCGainMatch::Init(TTree *tree){
 
 }
 
-Bool_t PCGainMatch::Notify(){
+Bool_t gainmatch::Notify(){
 
    return kTRUE;
 }
 
-void PCGainMatch::SlaveBegin(TTree * /*tree*/){
+void gainmatch::SlaveBegin(TTree * /*tree*/){
 
    TString option = GetOption();
 
 }
 
-void PCGainMatch::SlaveTerminate(){
+void gainmatch::SlaveTerminate(){
 
 }
 
 
-#endif // #ifdef Analyzer_cxx
+#endif // #ifdef gainmatch_cxx

mapping.h

@@ -19,7 +19,6 @@ const std::map<int, unsigned short> board = {
    {4, 22129},
    {5, 15529},
    {6, 15528},
-  //  {7,89},
    {7,   334},
    {8,   379},
    {9,   325},
@@ -28,14 +27,13 @@ const std::map<int, unsigned short> board = {
 const int nBd = board.size();
 
 const int nV1740 = 7;
-const int nV1725 = 4;
+const int nV1725 = 3;
 
 //+++++++++++++++++++ detID;
 // The detectors are seperated into 2 type: SuperX3, QQQ, and PC
 // the SuperX3 has 24 detectors for each kind, wach detector has 12 channels
 // the QQQ has 4 detectors for each kind, each detector has 32 channels
 // the PC has 2 types, anode and cathode, each has 24 channels
-// the MISC has 6 channels, the lollipop IC and siliscon followed by the hotneedle IC, as well as the Rf and MCP
 // The detID = Type * 10000 + index * 100 + channel
 // fro example, detID(superX3-8, ch-7) = 00807
 
@@ -78,23 +76,17 @@ const std::vector<int> mapping = {
     -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,
     -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,
     -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,
- //================== 89
-  //  0     1      2      3      4      5      6      7      8      9     10     11     12     13     14     15
-//  30004,    -1, 30003,    -1,    -1,    -1,    -1,    -1,     -1,   -1,    -1,    -1,    -1,    -1,    -1,    -1,
   //================== 334
- 20116, 20117, 20118, 20119,    -1, 20121, 20122, 20123, 20016, 20017, 20018,    -1, 20020, 20021, 20022, 20023,
+ 20116, 20117, 20118, 20119, 20120, 20121, 20122, 20123, 20016, 20017, 20018, 20019, 20020, 20021, 20022, 20023,
   //================== 379
-  -1  , 20001, 20002, 20003, 20004, 20005,    -1, 20007, 20008,    -1, 20010, 20011, 20012, 20013, 20014, 20015,
+ 20000, 20001, 20002, 20003, 20004, 20005,    -1, 20007, 20008,    -1, 20010, 20011, 20012, 20013, 20014, 20015,
   //================== 325
- 20100, 20101, 20102, 20103, 20104, 20105, 20106, 20107, 20108, 20109, 20110, 20111, 20112,    -1, 20114, 20115,
+ 20100, 20101, 20102, 20103, 20104, 20105, 20106, 20107, 20108, 20109, 20110, 20111, 20112, 20113, 20114, 20115,
   //================== 405
-  //  0     1      2      3      4      5      6      7      8      9     10     11     12     13     14     15
+ 20006,    -1,    -1, 20009,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1
- 20006,    -1, 30005, 20009,    -1, 20120, 20000, 20019, 20113, 30000, 30004, 30001, 30002,    -1, 30003,   -1
+
 };
 
-//MCP moved from channel 1 to 2 after Run number 322
-//MCP and Rf moved to ch 0 and 1 after Run number after Run282
-//moved back to ch 
 void PrintMapping(){
 
   int digiID = 0;
@@ -149,9 +141,8 @@ void PrintMapping(){
 
         printf("\033[35m%3d(%2d)\033[0m|", detID, ch);
         
-      }else if( typeID == 3){ // MISC
+      }else{
 
-        printf("\033[33m%3d(%2d)\033[0m|", detID, ch);
 
       }
     }
@@ -224,6 +215,7 @@ void GenMapping(std::string mapFile){
       detID += 20000;
       if( words[3] == "ANODE") detID += atoi(words[4].c_str());
       if( words[3] == "CATHODE") detID += 100 + atoi(words[4].c_str());
+
     }
 
     if( words[2] == "blank") {