modified: .gitignore

modified:   .vscode/settings.json
	modified:   Analyzer.C
	modified:   Calibration.C
	modified:   GainMatchQQQ.C
	modified:   QQQ_Calcheck.C

.gitignore

@@ -4,9 +4,23 @@ EventBuilder*
 *.pcm
 *.root
 *.exe
-
+*.txt
 Mapper
 AnasenMS
 
 data/
 data_proton/
+Sudarshan/
+Analyzer_C_ACLiC_dict0713aaa966_dictContent.h
+.gitignore
+Analyzer_C_ACLiC_dict5411fecd5c_dictUmbrella.h
+gainmatch.C
+gainmatch.h
+MakePlotsQQQ.C
+MakePlotsQQQ.h
+MakePlotsSX3.C
+MakePlotsSX3.h
+qqq_gains_det3.dat
+qqq_relative_gains.dat
+Armory/CorrelateQQQ.h
+QQQStage2.C

.vscode/c_cpp_properties.json

@@ -1,23 +1,11 @@
 {
     "configurations": [
         {
-            "name": "splitpole",
+            "name": "Hades",
             "includePath": [
                 "${workspaceFolder}/**",
-                "/home/splitpole/cern/root/**"
-            ],
-            "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"
+                "/usr/include/x86_64-linux-gnu/qt6/**",
+                "/usr/local/cern/root_v6.26.06/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,38 @@
             "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",
+            "cStandard": "c17",
+            "cppStandard": "gnu++17",
+            "intelliSenseMode": "linux-gcc-x64"
+        },
+        {
+            "name": "VigneshROG",
+            "includePath": [
+                "${workspaceFolder}/**",
+                "/home/vsitaraman/root/include/**"
             ],
             "defines": [],
             "compilerPath": "/usr/bin/gcc",

.vscode/settings.json

@@ -90,6 +90,41 @@
     "processRun.C": "cpp",
     "TrackRecon.C": "cpp",
     "processRuns.C": "cpp",
-    "Analysis.C": "cpp"
-  }
+    "Analysis.C": "cpp",
+    "datastructs.h": "c",
+    "ANASENPlotEdit.C": "cpp",
+    "GetMean_Q3_new.C": "cpp",
+    "AlphaCal_new.C": "cpp",
+    "f1.C": "cpp",
+    "GeoCal_Maria_new.C": "cpp",
+    "PCPulser_All_new.C": "cpp",
+    "PosCal_2.C": "cpp",
+    "AutoFit.C": "cpp",
+    "Fitting.C": "cpp",
+    "PCGainMatch.C": "cpp",
+    "Analyzer1.C": "cpp",
+    "FitHistogramsWithTSpectrum_Sequential_Improved.C": "cpp",
+    "PlotAndFitCentroids.C": "cpp",
+    "MatchAndPlotCentroids.C": "cpp",
+    "GainMatch.C": "cpp",
+    "GainMatchSX3.C": "cpp",
+    "RelBack_Fix_new.C": "cpp",
+    "SiRelativeGains_Step1_new.C": "cpp",
+    "charconv": "cpp",
+    "format": "cpp",
+    "GainMatchSX3Front.C": "cpp",
+    "GainMatchSX3Front1.C": "cpp",
+    "Calibration.C": "cpp",
+    "GainMatchQQQ.C": "cpp",
+    "UTF-8gainmatch.C": "cpp",
+    "MakePlotsQQQ.C": "cpp",
+    "MakePlotsSX3.C": "cpp",
+    "QQQ_Calibcheck.C": "cpp",
+    "QQQ_Calcheck.C": "cpp",
+    "makeplots.C": "cpp",
+    "GlobalMinimizeQQQ.C": "cpp",
+    "QQQStage2.C": "cpp",
+    "inspect.C": "cpp"
+  },
+  "github-enterprise.uri": "https://fsunuc.physics.fsu.edu"
 }

Analyzer.h

@@ -18,6 +18,7 @@ public :
    Det sx3;
    Det qqq;
    Det pc ;
+   Det misc;
 
    ULong64_t       evID;
    UInt_t          run;
@@ -40,6 +41,13 @@ 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() { }
@@ -92,6 +100,13 @@ 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

@@ -0,0 +1,402 @@
+#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");
+
+}

Analyzer1.h

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

Armory/ANASEN_model.C

@@ -103,8 +103,8 @@ void ANASEN_model(int anodeID1 = -1, int anodeID2 = -1, int cathodeID1 = -1, int
                                                        new TGeoRotation("rot1", 360/nSX3 * (i + 0.5), 0., 0.)));
   }
 
-  const int qqqR1 = 10;
-  const int qqqR2 = 50;
+  const int qqqR1 = 50;
+  const int qqqR2 = 100;
   TGeoVolume *qqq = geom->MakeTubs("qqq", Al, qqqR1, qqqR2, 0.5, 5, 85);
   qqq->SetLineColor(7);
   for( int i = 0; i < 4; i++){

Armory/ClassPW.h

@@ -4,15 +4,18 @@
 #include <cstdio>
 #include <TMath.h>
 #include <TVector3.h>
+#include <TRandom.h>
 
-struct PWHitInfo{
+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(){
+  void Clear()
+  {
     nearestWire.first = -1;
     nearestWire.second = -1;
     nearestDist.first = 999999999;
@@ -24,50 +27,65 @@ struct PWHitInfo{
   }
 };
 
-//!########################################################
-class PW{ // proportional wire
+struct Coord
+{
+  float x, y, z;
+  Coord() : x(0), y(0), z(0) {}
+  Coord(const TVector3 &vec)
+  {
+    x = vec.X(); // TVector3's X() returns the x-coordinate
+    y = vec.Y(); // TVector3's Y() returns the y-coordinate
+    z = vec.Z(); // TVector3's Z() returns the z-coordinate
+  }
+};
+
+//! ########################################################
+class PW
+{ // proportional wire
 public:
-  PW(){ ClearHitInfo();};
-  ~PW(){};
+  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;}
+  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();}
+  std::vector<std::pair<TVector3, TVector3>> An; // the anode wire position vector in space
+  std::vector<std::pair<TVector3, TVector3>> Ca; // the cathode wire position vector in space
+
+  TVector3 GetTrackPos() const { return trackPos; }
+  TVector3 GetTrackVec() const { return trackVec; }
+  double GetTrackTheta() const { return trackVec.Theta(); }
+  double GetTrackPhi() const { return trackVec.Phi(); }
   double GetZ0();
 
-  int GetNumWire() const {return nWire;}
-  double GetDeltaAngle() const {return dAngle;}
-  double GetAnodeLength() const {return anodeLength;}
-  double GetCathodeLength() const {return cathodeLength;}
-  TVector3 GetAnodeDn(short id) const {return An[id].first;}
-  TVector3 GetAnodeUp(short id) const {return An[id].second;}
-  TVector3 GetCathodeDn(short id) const {return Ca[id].first;}
-  TVector3 GetCathodeUp(short id) const {return Ca[id].second;}
+  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 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();}
+  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 CalTrack2(TVector3 sx3Pos, TVector3 anodeInt, bool verbose = false);
 
-  double CircularMean(std::vector<std::pair<int, double>> wireList);
-
-  void Print(){
+  void Print()
+  {
     printf("     The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first,
            hitInfo.nearestDist.first,
            hitInfo.nearestWire.second,
@@ -80,7 +98,6 @@ public:
   }
 
 private:
-
   PWHitInfo hitInfo;
 
   TVector3 trackPos;
@@ -88,7 +105,7 @@ private:
 
   const int nWire = 24;
   const int wireShift = 3;
-  const float zLen = 380; //mm
+  const float zLen = 380; // mm
   const float radiusA = 37;
   const float radiusC = 43;
 
@@ -96,23 +113,25 @@ private:
   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 
+  // 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){
+  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));
+    return TMath::Abs(nd.Dot(a1 - b2));
   }
-
 };
 
-inline void PW::ClearHitInfo(){
+inline void PW::ClearHitInfo()
+{
   hitInfo.Clear();
 }
 
-inline void PW::ConstructGeo(){
+inline void PW::ConstructGeo()
+{
 
   An.clear();
   Ca.clear();
@@ -120,108 +139,137 @@ inline void PW::ConstructGeo(){
   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.
+  // 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++ ){
+  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);
+    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);
+    // Cathod rotate left-hand with the 3 wire offset accounted for (+1 from the calculated offset from the PC coincidence spectrum)
+    q1.first.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()),
+                    radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()),
+                    zLen / 2);
+    q1.second.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i + 1) + TMath::PiOver2()),
+                     radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i + 1) + TMath::PiOver2()),
+                     -zLen / 2);
     Ca.push_back(q1);
   }
+  // correcting for the fact that the order of the cathode wires is reversed
+  std::reverse(Ca.begin(), Ca.end());
+  // adjusting for the 3 wire offset, the rbegin and rend are used as the rotation of the wires is done in the opposite direction i.e. 1,2,3 -> 3,1,2
+  // NOT NECESSARY ANY MORE, HAS BEEN IMCORPORATED INTO THE WIREOFFSET IN THE BEGINNING
+  // std::rotate(Ca.rbegin(), Ca.rbegin() + 4, Ca.rend());
 
   dAngle = wireShift * TMath::TwoPi() / nWire;
-  anodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusA * TMath::Sin(dAngle/2),2) );
-  cathodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusC * TMath::Sin(dAngle/2),2) );
+  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 ){
+inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
+{
 
   hitInfo.Clear();
   double phi = direction.Phi();
 
-  for( int i = 0; i < nWire; i++){
+  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 (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 ){
+    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();
+    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 (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 ){
+    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);
+    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 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){
+  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{
+  }
+  else
+  {
     hitInfo.nextNearestWire.first = aaa2;
     hitInfo.nextNearestDist.first = haha2;
   }
 
   short cathode1 = hitInfo.nearestWire.second;
-  short ccc1 = cathode1 - 1; if( ccc1 < 0 ) ccc1 += nWire;
+  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){
+  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{
+  }
+  else
+  {
     hitInfo.nextNearestWire.second = ccc2;
     hitInfo.nextNearestDist.second = haha2;
   }
 
-  if( verbose ) Print();
+  if (verbose)
+    Print();
 }
 
-inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose){
+inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose)
+{
 
   trackPos = sx3Pos;
 
@@ -231,72 +279,37 @@ inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbo
   // 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()); 
-
+  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;
+inline void PW::CalTrack2(TVector3 siPos, TVector3 anodeInt, bool verbose)
+{
 
-  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()); 
+  float mx, my;
+  double z;
+  mx = siPos.X() / (siPos.X() - anodeInt.X());
+  my = siPos.Y() / (siPos.Y() - anodeInt.Y());
+  z=siPos.Z() + mx * (anodeInt.Z() - siPos.Z());
+  // if (mx == my)
+  {
+    trackVec=TVector3(0,0,z);
+  }
 
+  if (verbose)
+    printf("X slope = %f and Y slope = %f \n", mx, my);
 }
 
-inline double PW::GetZ0(){
+inline double PW::GetZ0()
+{
 
   double x = trackPos.X();
   double y = trackPos.Y();
-  double rho = TMath::Sqrt(x*x + y*y);
+  double rho = TMath::Sqrt(x * x + y * y);
   double theta = trackVec.Theta();
 
-  return trackPos.Z() - rho / TMath::Tan(theta);
-
-}
-
-inline double PW::CircularMean(std::vector<std::pair<int, double>> wireList){
-
-  //use unit vector, wireID start from Zero
-  double xCom = 0, yCom = 0;
-  for( size_t i = 0; i < wireList.size() ; i++){
-    xCom += TMath::Cos(TMath::TwoPi() * wireList[i].first / nWire) * wireList[i].second;
-    yCom += TMath::Sin(TMath::TwoPi() * wireList[i].first / nWire) * wireList[i].second;
-  }
-
-  //calculate the angle of the summed unit vectors
-  double angle = TMath::ATan2(yCom, xCom);
-  if( angle < 0 ) angle += TMath::TwoPi(); // convert the angle from 0 to 2 pi
-
-  return angle/ TMath::TwoPi() * nWire;
-
+  return trackVec.Z();
 }
 
 #endif

Armory/HistPlotter.h

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

Armory/LICENSE

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

Armory/README.md

@@ -0,0 +1,6 @@
+# HistPlotter
+ - Header-only class to encapsulate CERN ROOT 1D/2D Histogram plotting and application of TCuts.
+ - Can specify folder hierarchy while setting up fills, currently supports only one level.
+ - Cuts specified using a two-column text file containing cut names, and target .root files. These .root files must contain a TCutG of name "CUTG".
+ - Tested for use in macros, with TSelector design pattern and compiled code.
+ - Idea inspired from MyFill() pattern created by github user gwm17

Calibration.C

@@ -0,0 +1,270 @@
+
+#define Calibration_cxx
+
+#include <TH2.h>
+#include <TF1.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TCutG.h>
+#include <fstream>
+#include <utility>
+#include <algorithm>
+#include "Armory/HistPlotter.h"
+#include "TVector3.h"
+#include "Calibration.h"
+
+TH2F *hQQQFVB;
+HistPlotter *plotter;
+int padID = 0;
+
+TCutG *cut;
+std::map<std::tuple<int, int, int>, std::vector<std::pair<double, double>>> dataPoints;
+
+bool qqqEcut = false;
+
+// Gain Arrays
+const int MAX_QQQ = 4;
+const int MAX_RING = 16;
+const int MAX_WEDGE = 16;
+double qqqwGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
+// double qqqrGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
+bool qqqwGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
+// bool qqqrGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
+
+void Calibration::Begin(TTree * /*tree*/)
+{
+    plotter = new HistPlotter("Calib.root", "TFILE");
+    // ----------------------- Load QQQ Gains
+    {
+        std::string filename = "qqq_GainMatch.txt";
+        std::ifstream infile(filename);
+        if (!infile.is_open())
+        {
+            std::cerr << "Error opening " << filename << "!" << std::endl;
+        }
+        else
+        {
+            int det, ring, wedge;
+            double gainw, gainr;
+            while (infile >> det >> ring >> wedge >> gainw >> gainr)
+            {
+                qqqwGain[det][ring][wedge] = gainw;
+                // qqqrGain[det][ring][wedge] = gainr;
+                qqqwGainValid[det][ring][wedge] = (gainw > 0);
+                // qqqrGainValid[det][ring][wedge] = (gainr > 0);
+            }
+            infile.close();
+            std::cout << "Loaded QQQ gains from " << filename << std::endl;
+        }
+    }
+    for (int det = 0; det < MAX_QQQ; det++)
+    {
+        for (int ring = 0; ring < MAX_RING; ring++)
+        {
+            for (int wedge = 0; wedge < MAX_WEDGE; wedge++)
+            {
+                TString hname = Form("hCal_qqq%d_ring%d_wedge%d", det, ring, wedge);
+                TString htitle = Form("QQQ det%d ring%d wedge%d; Energy (arb); Counts", det, ring, wedge);
+                // hQQQSpectra[det][ring][wedge] = new TH1F(hname, htitle, 4000, 0, 16000);
+            }
+        }
+    }
+}
+
+Bool_t Calibration::Process(Long64_t entry)
+{
+    b_qqqMulti->GetEntry(entry);
+    b_qqqID->GetEntry(entry);
+    b_qqqCh->GetEntry(entry);
+    b_qqqE->GetEntry(entry);
+    b_qqqT->GetEntry(entry);
+
+    qqq.CalIndex();
+
+    for (int i = 0; i < qqq.multi; i++)
+    {
+        for (int j = i + 1; j < qqq.multi; j++)
+        {
+            if (qqq.e[i] > 100)
+                qqqEcut = true;
+            if (qqq.id[i] == qqq.id[j])
+            {
+                int chWedge = -1;
+                int chRing = -1;
+                float eWedgeRaw = 0.0;
+                float eWedge = 0.0;
+                float eRingRaw = 0.0;
+                float eRing = 0.0;
+                if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && /*qqqrGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16] &&*/ qqqwGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
+                {
+                    chWedge = qqq.ch[i];
+                    eWedgeRaw = qqq.e[i];
+                    eWedge = qqq.e[i] * qqqwGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
+                    // printf("Wedge E: %.2f  Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]);
+                    chRing = qqq.ch[j] - 16;
+                    eRingRaw = qqq.e[j];
+                    eRing = qqq.e[j];// * qqqrGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
+                }
+                else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && /*qqqrGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16] &&*/ qqqwGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
+                {
+                    chWedge = qqq.ch[j];
+                    eWedge = qqq.e[j] * qqqwGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
+                    eWedgeRaw = qqq.e[j];
+
+                    chRing = qqq.ch[i] - 16;
+                    eRing = qqq.e[i];// * qqqrGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
+                    eRingRaw = qqq.e[i];
+                }
+                else
+                    continue;
+
+                // hQQQFVB->Fill(eWedge, eRing);
+                plotter->Fill2D(Form("hRaw_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000, eWedgeRaw, eRingRaw, "ERaw");
+                plotter->Fill2D(Form("hGM_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000, eWedge, eRing, "EGM");
+                plotter->Fill2D("hRawQQQ", 4000, 0, 16000, 4000, 0, 16000, eWedgeRaw, eRingRaw);
+                plotter->Fill2D("hGMQQQ", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing);
+
+                TString histName = Form("hQQQFVB_id%d_r%d_w%d", qqq.id[i], chRing, chWedge);
+                // TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
+                // if (!hist2d)
+                // {
+                //     hist2d = new TH2F(histName, Form("QQQ Det%d R%d W%d;Wedge E;Ring E", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000);
+                // }
+
+                // hist2d->Fill(eWedge, eRing);
+                // if (cut && cut->IsInside(eWedge, eRing))
+                const double MIN_ADC = 1500.0;
+                const double MAX_ADC = 3000.0;
+
+                // if (eWedge >= MIN_ADC && eWedge <= MAX_ADC &&
+                //     eRing >= MIN_ADC && eRing <= MAX_ADC)
+                double ratio = (eWedge > 0.0) ? (eRing / eWedge) : 0.0;
+
+                double maxslope = 1.5;
+
+                bool validPoint = false;
+                if (ratio < maxslope && ratio > 1. / maxslope)
+                {
+                    // Accumulate data for gain matching
+                    dataPoints[{qqq.id[i], chRing, chWedge}].emplace_back(eWedge, eRing);
+                }
+            }
+        }
+    }
+
+    return kTRUE;
+}
+
+void Calibration::Terminate()
+{
+    const double AM241_PEAK = 5485.56;
+    const double P_PEAK = 7000; // keV
+
+    double calibArray[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
+    bool calibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
+
+    std::ofstream outFile("qqq_Calib.txt");
+    if (!outFile.is_open())
+    {
+        std::cerr << "Error opening qqq_Calib.txt!" << std::endl;
+        return;
+    }
+
+    //----------------------------------------------------------------------
+    // 1. Create per–channel 1D spectra in ADC from stored gain-matched data
+    //----------------------------------------------------------------------
+
+    std::map<std::tuple<int, int, int>, TH1F *> spectra;
+
+    for (auto &kv : dataPoints)
+    {
+        int det, ring, wedge;
+        std::tie(det, ring, wedge) = kv.first;
+
+        TString hname = Form("hSpec_d%d_r%d_w%d", det, ring, wedge);
+        TH1F *h = new TH1F(hname, hname, 4000, 0, 16000);
+
+        for (auto &p : kv.second)
+        {
+            double eWedge = p.first; // already gain-matched ADC
+            double eRing = p.second;
+
+            // Use ring ADC for calibration (cleaner alpha peak)
+            h->Fill(eRing);
+        }
+
+        spectra[kv.first] = h;
+    }
+
+    //----------------------------------------------------------------------
+    // 2. Fit Am-241 peak and extract keV/ADC calibration slope
+    //----------------------------------------------------------------------
+
+    for (auto &kv : spectra)
+    {
+        int det, ring, wedge;
+        std::tie(det, ring, wedge) = kv.first;
+        TH1F *h = kv.second;
+
+        if (!h || h->GetEntries() < 50)
+            continue;
+
+        int binMax = h->GetMaximumBin();
+        double adcPeak = h->GetXaxis()->GetBinCenter(binMax);
+
+        if (adcPeak <= 0)
+            continue;
+
+        // double slope_keV = AM241_PEAK / adcPeak; // keV per ADC
+        double slope_keV = P_PEAK / adcPeak; // keV per ADC
+
+        calibArray[det][ring][wedge] = slope_keV;
+        calibValid[det][ring][wedge] = true;
+
+        outFile << det << " " << ring << " " << wedge << " "
+                << slope_keV << "\n";
+
+        // printf("QQQ DET=%d R=%d W=%d  ADCpeak=%.1f  slope_keV=%.6f\n",det, ring, wedge, adcPeak, slope_keV);
+    }
+
+    outFile.close();
+    std::cout << "Wrote QQQ calibration file qqq_Calib.txt\n";
+
+    //----------------------------------------------------------------------
+    // 3. Build fully calibrated 2D combined histogram
+    //----------------------------------------------------------------------
+
+    TH2F *hCal = new TH2F("hCal",
+                          "All QQQ Calibrated;Wedge Energy (keV);Ring Energy (keV)",
+                          800, 0, 7000,
+                          800, 0, 7000);
+
+    for (auto &kv : dataPoints)
+    {
+        int det, ring, wedge;
+        std::tie(det, ring, wedge) = kv.first;
+
+        if (!calibValid[det][ring][wedge])
+            continue;
+
+        double slope = calibArray[det][ring][wedge];
+
+        for (auto &p : kv.second)
+        {
+            double eWGM = p.first; // gain matched ADC
+            double eRGM = p.second;
+
+            double eWkeV = eWGM * slope / 1000;
+            double eRkeV = eRGM * slope / 1000;
+
+            hCal->Fill(eWkeV, eRkeV);
+            plotter->Fill2D("hCalQQQ", 4000, 0, 10, 4000, 0, 10, eWkeV, eRkeV);
+            plotter->Fill2D(Form("hRCal_qqq%d", det), 16, 0, 15, 400, 0, 24, ring, eRkeV, "RingCal");
+            plotter->Fill2D(Form("hWCal_qqq%d", det), 16, 0, 15, 400, 0, 24, wedge, eWkeV, "WedgeCal");
+        }
+    }
+
+    plotter->FlushToDisk();
+    std::cout << "Calibrated 2D QQQ histogram saved.\n";
+}

Calibration.h

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

FitHistogramsWithTSpectrum_Sequential_Improved.C

@@ -0,0 +1,124 @@
+#include <TFile.h>
+#include <TH1.h>
+#include <TSpectrum.h>
+#include <TF1.h>
+#include <TCanvas.h>
+#include <vector>
+#include <iostream>
+#include <algorithm>
+#include <fstream>
+#include <TText.h>
+
+void FitHistogramsWithTSpectrum_Sequential_Improved() {
+    TFile *inputFile = new TFile("Histograms_anodes.root", "READ");
+    if (!inputFile || inputFile->IsZombie()) {
+        std::cerr << "Error opening the input file!" << std::endl;
+        return;
+    }
+
+    TCanvas *c1 = new TCanvas("c1", "Histogram Viewer", 800, 600);
+
+    // Open the output ASCII file to save the centroids
+    std::ofstream outFile("centroids.txt");
+    if (!outFile.is_open()) {
+        std::cerr << "Error opening output file!" << std::endl;
+        return;
+    }
+    outFile << "HistogramIndex\tPeakNumber\tCentroid\tAmplitude\tSigma" << std::endl;
+
+    for (int i = 0; i < 24; ++i) {
+        TH1 *histogram = dynamic_cast<TH1*>(inputFile->Get(Form("hCathode_%d", i)));
+        if (!histogram) {
+            std::cerr << "Failed to retrieve histogram_" << i << " from the file." << std::endl;
+            continue;
+        }
+
+        // Set range for peak search
+        double minX = 700;
+        double maxX = 25000;
+        histogram->GetXaxis()->SetRangeUser(minX, maxX);
+
+        // Draw the histogram
+        c1->cd();
+        histogram->Draw();
+
+        // Peak search using TSpectrum
+        const int maxPeaks = 5;
+        TSpectrum spectrumFinder(maxPeaks);
+        int nFound = spectrumFinder.Search(histogram, 2, "", 0.01);
+
+        if (nFound <= 0) {
+            std::cerr << "No peaks found for histogram " << i << std::endl;
+            continue;
+        }
+
+        Double_t *xPositions = spectrumFinder.GetPositionX();
+        Double_t *yPositions = spectrumFinder.GetPositionY();
+        std::vector<std::pair<Double_t, Double_t>> peaks;
+
+        // Collect and sort peaks by X position
+        for (int j = 0; j < nFound; ++j) {
+            peaks.emplace_back(xPositions[j], yPositions[j]);
+        }
+        std::sort(peaks.begin(), peaks.end());
+
+        // Fit each peak with a Gaussian
+        for (int j = 0; j < peaks.size(); ++j) {
+            Double_t peakX = peaks[j].first;
+            Double_t peakY = peaks[j].second;
+            Double_t initialAmplitude = peakY;  // Better initial guess
+            Double_t initialCentroid = peakX;   // Centroid based on peak position
+            Double_t initialSigma = 60.0;
+            // Define Gaussian with initial parameters
+            TF1 *gaussFit = new TF1(Form("gauss_%d", j), "gaus", peakX - 200, peakX + 200);
+            //gaussFit->SetParameters(peakY, peakX, 25.0); // Initial guesses for amplitude, mean, sigma
+            gaussFit->SetParameters(initialAmplitude, initialCentroid, initialSigma);
+            // Perform fit
+            int fitStatus = histogram->Fit(gaussFit, "RQ+");
+            if (fitStatus != 0) {
+                std::cerr << "Fit failed for peak " << j + 1 << " in histogram " << i << std::endl;
+                delete gaussFit;
+                continue;
+            }
+
+            // Retrieve fit parameters
+            double amplitude = gaussFit->GetParameter(0);
+            double centroid = gaussFit->GetParameter(1);
+            double sigma = gaussFit->GetParameter(2);
+            double amplitudeError = gaussFit->GetParError(0);
+            double centroidError = gaussFit->GetParError(1);
+            double sigmaError = gaussFit->GetParError(2);
+
+            // Chi-squared value
+            double chi2 = gaussFit->GetChisquare();
+            int ndf = gaussFit->GetNDF();
+            outFile << i << "\t" << j + 1 << "\t" << centroid << std::endl;
+            gaussFit->SetLineColor(kRed);
+            gaussFit->Draw("SAME");
+            TText *text = new TText();
+            text->SetNDC();
+            text->SetTextSize(0.03);
+            text->SetTextColor(kRed);
+            //text->DrawText(0.15, 0.8 - j * 0.05, Form("Peak %d: Amp=%.2f, Mean=%.2f, Sigma=%.2f", j + 1, amplitude, centroid, sigma));
+            text->DrawText(0.15, 0.8 - j * 0.05, 
+                Form("Peak %d: Amp=%.2f±%.2f, Mean=%.2f±%.2f, Sigma=%.2f±%.2f, Chi2/NDF=%.2f",
+                     j + 1, amplitude, amplitudeError, centroid, centroidError, sigma, sigmaError, chi2 / ndf));
+            // Save results
+            
+            
+            // Clean up
+            delete gaussFit;
+        }
+
+        // Update canvas for visualization
+        c1->Update();
+        std::cout << "Press Enter to view the next histogram..." << std::endl;
+        c1->WaitPrimitive(); // Wait until Enter is pressed in the ROOT console
+    }
+
+    // Close resources
+    inputFile->Close();
+    outFile.close();
+    delete c1;
+}
+

GainMatchQQQ.C

@@ -0,0 +1,309 @@
+#define GainMatchQQQ_cxx
+
+#include "GainMatchQQQ.h"
+#include <TH2.h>
+#include <TF1.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TCutG.h>
+#include <fstream>
+#include <utility>
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+#include "Armory/HistPlotter.h"
+#include "TVector3.h"
+#include "TGraphErrors.h"
+#include "TF1.h"
+#include <cmath>
+
+TH2F *hQQQFVB;
+HistPlotter *plotter;
+
+int padID = 0;
+
+TCutG *cut;
+std::map<std::tuple<int, int, int>, std::vector<std::pair<double, double>>> dataPoints;
+
+void GainMatchQQQ::Begin(TTree * /*tree*/)
+{
+    plotter = new HistPlotter("GainQQQ.root", "TFILE");
+    TString option = GetOption();
+
+    hQQQFVB = new TH2F("hQQQFVB", "QQQ Front vs Back; Front E; Back E", 800, 0, 16000, 800, 0, 16000);
+
+    // Load the TCutG object
+    TFile *cutFile = TFile::Open("qqqcorr.root");
+    if (!cutFile || cutFile->IsZombie())
+    {
+        std::cerr << "Error: Could not open qqqcorr.root" << std::endl;
+        return;
+    }
+    cut = dynamic_cast<TCutG *>(cutFile->Get("qqqcorr"));
+    if (!cut)
+    {
+        std::cerr << "Error: Could not find TCutG named 'qqqcorr' in qqqcorr.root" << std::endl;
+        return;
+    }
+    cut->SetName("qqqcorr"); // Ensure the cut has the correct name
+}
+
+Bool_t GainMatchQQQ::Process(Long64_t entry)
+{
+
+    int ringMults[16] = {0};
+    int wedgeMults[16] = {0};
+    std::vector<std::tuple<int, int, int, double, double>> events;
+
+    b_qqqMulti->GetEntry(entry);
+    b_qqqID->GetEntry(entry);
+    b_qqqCh->GetEntry(entry);
+    b_qqqE->GetEntry(entry);
+    b_qqqT->GetEntry(entry);
+
+    qqq.CalIndex();
+
+    for (int i = 0; i < qqq.multi; i++)
+    {
+        for (int j = i + 1; j < qqq.multi; j++)
+        {
+            if (qqq.id[i] == qqq.id[j])
+            {
+                int chWedge = -1;
+                int chRing = -1;
+                float eWedge = 0.0;
+                float eRing = 0.0;
+                if (qqq.ch[i] < 16 && qqq.ch[j] >= 16)
+                {
+                    chWedge = qqq.ch[i];
+                    eWedge = qqq.e[i];
+                    chRing = qqq.ch[j] - 16;
+                    eRing = qqq.e[j];
+                }
+                else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16)
+                {
+                    chWedge = qqq.ch[j];
+                    eWedge = qqq.e[j];
+                    chRing = qqq.ch[i] - 16;
+                    eRing = qqq.e[i];
+                }
+                else
+                    continue;
+                ringMults[chRing]++;
+                wedgeMults[chWedge]++;
+                hQQQFVB->Fill(eWedge, eRing);
+                events.emplace_back(qqq.id[i], chRing, chWedge, eRing, eWedge);
+                plotter->Fill2D(Form("hRaw_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 800, 0, 3000, 800, 0, 3000, eWedge, eRing, "hRawQQQ");
+                // double ratio = (eWedge > 0.0) ? (eRing / eWedge) : 0.0;
+                // double maxslope = 1.5;
+
+                // bool validPoint = false;
+                // if (ratio < maxslope && ratio > 1. / maxslope)
+                // {
+                //     // Accumulate data for gain matching
+                //     dataPoints[{qqq.id[i], chRing, chWedge}].emplace_back(eWedge, eRing);
+                //     plotter->Fill2D("hAll_in", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing);
+                //     validPoint = true;
+                // }
+
+                // if (!validPoint)
+                // {
+                //     plotter->Fill2D("hAll_out", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing);
+                // }
+            }
+        }
+    }
+
+    for (auto tuple : events)
+    {
+        auto [id, chr, chw, er, ew] = tuple;
+        if (ringMults[chr] > 1 || wedgeMults[chw] > 1)
+            continue; // ignore multiplicity > 1 events
+        double ratio = (ew > 0.0) ? (er / ew) : 0.0;
+        double maxslope = 1.5;
+
+        bool validPoint = false;
+        if (ratio < maxslope && ratio > 1. / maxslope)
+        {
+                // Accumulate data for gain matching
+                dataPoints[{id, chr, chw}].emplace_back(ew, er);
+                plotter->Fill2D("hAll_in", 4000, 0, 16000, 4000, 0, 16000, ew, er);
+                validPoint = true;
+                    }
+        if (!validPoint)
+        {
+            plotter->Fill2D("hAll_out", 4000, 0, 16000, 4000, 0, 16000, ew, er);
+        }
+    }
+
+    return kTRUE;
+}
+
+
+void GainMatchQQQ::Terminate()
+{
+    const int MAX_DET = 4;
+    const int MAX_RING = 16;
+    const int MAX_WEDGE = 16;
+
+    // We store gains locally just for the "corrected" plot, 
+    // but the file will output Slopes for the global minimizer.
+    double gainW[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}};
+    double gainR[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}};
+    bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}};
+
+    // Output file for the Minimizer
+    std::ofstream outFile("qqq_GainMatch.txt");
+    
+    // Benchmark/Debug file
+    std::ofstream benchFile("benchmark_diff.txt"); 
+    benchFile << "ID Wedge Ring Chi2NDF Slope SlopeErr" << std::endl;
+
+    if (!outFile.is_open()) { std::cerr << "Error opening output file!" << std::endl; return; }
+
+    const int MIN_POINTS = 50;
+    const int MAX_ITER = 3;        // Outlier rejection passes
+    const double CLIP_SIGMA = 2.5; // Sigma threshold for outliers
+
+    for (const auto &kv : dataPoints)
+    {
+        auto key = kv.first;
+        auto [id, ring, wedge] = key;
+        const auto &pts = kv.second;
+
+        if (pts.size() < (size_t)MIN_POINTS) continue;
+
+        std::vector<std::pair<double, double>> current_pts = pts;
+        
+        double finalSlope = 0.0;
+        double finalSlopeErr = 0.0;
+        bool converged = false;
+
+        // --- Iterative Fitting ---
+        for (int iter = 0; iter < MAX_ITER; ++iter)
+        {
+            if (current_pts.size() < (size_t)MIN_POINTS) break;
+
+            std::vector<double> x, y, ex, ey;
+
+            for (const auto &p : current_pts)
+            {
+                x.push_back(p.first);  // Wedge E
+                y.push_back(p.second); // Ring E
+                ex.push_back(std::sqrt(std::abs(p.first))); // Error in X (Poisson)
+                ey.push_back(std::sqrt(std::abs(p.second))); // Error in Y (Poisson)
+                
+                // Sanity check to avoid 0 error
+                if(ex.back() < 1.0) ex.back() = 1.0;
+                if(ey.back() < 1.0) ey.back() = 1.0;
+            }
+
+            // 2. Create Graph
+            TGraphErrors *gr = new TGraphErrors(current_pts.size(), x.data(), y.data(), ex.data(), ey.data());
+            
+            // 3. Fit Linear Function through Origin
+            TF1 *f1= new TF1("calibFit", "[0]*x", 0, 16000);
+            f1->SetParameter(0, 1.0);   
+
+            // "Q"=Quiet, "N"=NoDraw, "S"=ResultPtr
+            // We do NOT use "W" (Ignore weights), we want to use the errors we set.
+            int fitStatus = gr->Fit(f1, "QNS");
+
+            if (fitStatus != 0) {
+                delete gr; delete f1;
+                break; 
+            }
+
+            finalSlope = f1->GetParameter(0);
+            double chi2 = f1->GetChisquare();
+            double ndf = f1->GetNDF();
+            
+            // Get the statistical error on the slope
+            double rawErr = f1->GetParError(0);
+            
+            // SCALING ERROR:
+            // If Chi2/NDF > 1, the data scatters more than Poisson stats predict.
+            // // We inflate the error by sqrt(Chi2/NDF) to be conservative for the Minimizer.
+            // double redChi2 = (ndf > 0) ? (chi2 / ndf) : 1.0;
+            // double inflation = (redChi2 > 1.0) ? std::sqrt(redChi2) : 1.0;
+            
+            // finalSlopeErr = rawErr * inflation;
+
+            // 4. Outlier Rejection
+            if (iter == MAX_ITER - 1) {
+                converged = true;
+                delete gr; delete f1;
+                break;
+            }
+
+            // Calculate Residuals
+            std::vector<double> residuals;
+            double sumSqResid = 0.0;
+            for(size_t k=0; k<current_pts.size(); ++k) {
+                double val = f1->Eval(current_pts[k].first);
+                double res = current_pts[k].second - val;
+                residuals.push_back(res);
+                sumSqResid += res*res;
+            }
+        //     double sigma = std::sqrt(sumSqResid / current_pts.size());
+            
+        //     // Filter
+        //     std::vector<std::pair<double, double>> next_pts;
+        //     for(size_t k=0; k<current_pts.size(); ++k) {
+        //         if(std::abs(residuals[k]) < CLIP_SIGMA * sigma) {
+        //             next_pts.push_back(current_pts[k]);
+        //         }
+        //     }
+
+        //     if (next_pts.size() == current_pts.size()) {
+        //         converged = true; 
+        //         delete gr; delete f1;
+        //         break;
+        //     }
+        //     current_pts = next_pts;
+        //     delete gr; delete f1;
+        }
+
+        if (!converged || finalSlope <= 0) continue;
+
+        // --- Store/Output ---
+        
+        // 1. Save locally for the verification plot (hAll)
+        // Approximate local gain for plotting purposes only
+        double gW_local = std::sqrt(finalSlope);
+        double gR_local = 1.0 / gW_local;
+        gainW[id][ring][wedge] = gW_local;
+        gainR[id][ring][wedge] = gR_local;
+        gainValid[id][ring][wedge] = true;
+
+        // 2. Write to File for Minimizer
+        // Format: ID Wedge Ring Slope Error
+        outFile << id << " " << wedge << " " << ring << " " << finalSlope << " " << finalSlopeErr << std::endl;
+
+        // 3. Benchmark Info
+        benchFile << id << " " << wedge << " " << ring << " " 
+                  << finalSlope << " " << finalSlopeErr << std::endl;
+    }
+
+    outFile.close();
+    benchFile.close();
+    std::cout << "Gain matching with Errors complete." << std::endl;
+
+    // Plotting the corrected data (Visual check using local approx gains)
+    for (auto &kv : dataPoints)
+    {
+        int id, ring, wedge;
+        std::tie(id, ring, wedge) = kv.first;
+        if (!gainValid[id][ring][wedge]) continue;
+        auto &pts = kv.second;
+        for (auto &pr : pts)
+        {
+            double corrWedge = pr.first * gainW[id][ring][wedge];
+            double corrRing = pr.second * gainR[id][ring][wedge];
+            plotter->Fill2D("hAll", 4000, 0, 16000, 4000, 0, 16000, corrWedge, corrRing);
+        }
+    }
+
+    plotter->FlushToDisk();
+}

GainMatchQQQ.h

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

GainMatchSX3.C

@@ -0,0 +1,432 @@
+#define GainMatchSX3_cxx
+
+#include "GainMatchSX3.h"
+#include <TH2.h>
+#include <TF1.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TCutG.h>
+#include <fstream>
+#include <utility>
+#include <algorithm>
+#include <TProfile.h>
+#include "Armory/ClassSX3.h"
+#include "Armory/HistPlotter.h"
+#include <TGraphErrors.h>
+#include "TVector3.h"
+
+TH2F *hSX3FvsB;
+TH2F *hSX3FvsB_g;
+TH2F *hsx3IndexVE;
+TH2F *hsx3IndexVE_g;
+TH2F *hSX3;
+TH2F *hsx3Coin;
+
+int padID = 0;
+
+SX3 sx3_contr;
+TCutG *cut;
+TCutG *cut1;
+std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
+std::map<std::tuple<int, int, int, int>, int> comboCounts;
+
+const int MAX_DET = 24;
+const int MAX_UP = 4;
+const int MAX_DOWN = 4;
+const int MAX_BK = 4;
+
+double frontGainUp[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
+double frontGainDown[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
+bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
+TCanvas c("canvas", "canvas", 800, 600);
+
+// ==== Configuration Flags ====
+const bool interactiveMode = true; // If true: show canvas + wait for user
+const bool verboseFit = true;      // If true: print fit summary and chi²
+const bool drawCanvases = true;    // If false: canvases won't be drawn at all
+
+// HistPlotter plotter("SX3GainMatchBack.root");
+
+void GainMatchSX3::Begin(TTree * /*tree*/)
+{
+    TString option = GetOption();
+
+    hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
+    hSX3FvsB_g = new TH2F("hSX3FvsB_g", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
+    hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
+    hsx3IndexVE_g = new TH2F("hsx3IndexVE_g", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
+    hSX3 = new TH2F("hSX3", "SX3 Front v Back; Fronts; Backs", 8, 0, 8, 4, 0, 4);
+
+    hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12);
+
+    sx3_contr.ConstructGeo();
+
+    // Load the TCutG object
+    TFile *cutFile = TFile::Open("sx3cut.root");
+    if (!cutFile || cutFile->IsZombie())
+    {
+        std::cerr << "Error: Could not open sx3cut.root" << std::endl;
+        return;
+    }
+    cut = dynamic_cast<TCutG *>(cutFile->Get("sx3cut"));
+    if (!cut)
+    {
+        std::cerr << "Error: Could not find TCutG named 'sx3cut' in sx3cut.root" << std::endl;
+        return;
+    }
+    cut->SetName("sx3cut"); // Ensure the cut has the correct name
+
+    // Load the TCutG object
+    TFile *cutFile1 = TFile::Open("UvD.root");
+    bool cut1Loaded = (cut1 != nullptr);
+    cut1 = dynamic_cast<TCutG *>(cutFile1->Get("UvD"));
+    if (!cut1)
+    {
+        std::cerr << "Error: Could not find TCutG named 'UvD' in UvD.root" << std::endl;
+        return;
+    }
+    cut1->SetName("UvD");
+
+    // plotter.ReadCuts("cuts.txt");
+
+    std::string filename = "sx3_GainMatchfront.txt";
+    // std::string filename = "sx3_GainMatchfront.txt";
+
+    std::ifstream infile(filename);
+    if (!infile.is_open())
+    {
+        std::cerr << "Error opening " << filename << "!" << std::endl;
+        return;
+    }
+
+    int id, bk, u, d;
+    double gainup, gaindown;
+    while (infile >> id >> bk >> u >> d >> gainup >> gaindown)
+    {
+        frontGainUp[id][bk][u][d] = gainup;
+        frontGainDown[id][bk][u][d] = gaindown;
+        frontGainValid[id][bk][u][d] = true;
+        if(frontGainValid[id][bk][u][d]) {
+            // std::cout << "Loaded front gain for Det" << id << " Bk" << bk << " U" << u << " D" << d
+            //           << ": Up=" << gainup << ", Down=" << gaindown << std::endl;
+        }
+        else {
+            std::cout << "No valid front gain for Det" << id << " Bk" << bk << " U" << u << " D" << d << std::endl;
+        }
+    }
+}
+
+Bool_t GainMatchSX3::Process(Long64_t 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();
+
+    std::vector<std::pair<int, int>> ID;
+    for (int i = 0; i < sx3.multi; i++)
+    {
+
+        // for (int j = i + 1; j < sx3.multi; j++)
+        // {
+        //     if (sx3.id[i] == 3)
+        //         hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
+        // }
+        if (sx3.e[i] > 100)
+        {
+            ID.push_back(std::pair<int, int>(sx3.id[i], i));
+            hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
+        }
+    }
+
+    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; });
+
+        // start with the first entry in the sorted array: channels that belong to the same detector are together in sequenmce
+        std::vector<std::pair<int, int>> sx3ID;
+        sx3ID.push_back(ID[0]);
+        bool found = false;
+
+        for (size_t i = 1; i < ID.size(); i++)
+        { // Check if id of i belongs to the same detector and then add it to the detector ID vector
+            if (ID[i].first == sx3ID.back().first)
+            { // count the nunmber of hits that belong to the same detector
+                sx3ID.push_back(ID[i]);
+
+                if (sx3ID.size() >= 3)
+                {
+                    found = true;
+                }
+            }
+            else
+            { // the next event does not belong to the same detector, abandon the first event and continue with the next one
+                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;
+
+            // Build the correlated set once
+            for (size_t i = 0; i < sx3ID.size(); i++)
+            {
+                if (sx3.e[i] > 100)
+                {
+                    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] - 8;
+                        sx3EBk = sx3.e[index];
+                    }
+                }
+                
+            }
+            // Only if we found all three channels do we proceed
+            if (sx3ChUp >= 0 && sx3ChDn >= 0 && sx3ChBk >= 0)
+            {
+                // Fill once per correlated set
+                hSX3->Fill(sx3ChDn + 4, sx3ChBk);
+                hSX3->Fill(sx3ChUp, sx3ChBk);
+                hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
+                
+                if (frontGainValid[sx3ID[0].first][sx3ChBk][sx3ChUp / 2][sx3ChDn / 2])
+                {
+                    sx3EUp *= frontGainUp[sx3ID[0].first][sx3ChBk][sx3ChUp / 2][sx3ChDn / 2];
+                    sx3EDn *= frontGainDown[sx3ID[0].first][sx3ChBk][sx3ChUp / 2][sx3ChDn / 2];
+                }
+                else
+                {
+                    // printf("No front gain for Det%d Bk%d U%d D%d\n", sx3ID[0].first, sx3ChBk, sx3ChUp / 2, sx3ChDn / 2);
+                    sx3EUp = sx3EDn = 0.;
+                }
+                // plotter.Fill2D("hSX3F", 400, 0, 16000, 400, 0, 16000, sx3EUp + sx3EDn, sx3EBk);
+
+                // Pick detector ID from one of the correlated hits (all same detector)
+                int detID = sx3ID[0].first;
+
+                TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", detID, sx3ChUp, sx3ChDn, sx3ChBk);
+                TString histName1 = Form("UnCorr_id%d_U%d-D%dvsB%d", detID, sx3ChUp, sx3ChDn, sx3ChBk);
+
+                TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
+                TH2F *hist2d1 = (TH2F *)gDirectory->Get(histName1);
+                if (!hist2d)
+                {
+                    hist2d = new TH2F(histName, histName,
+                                      400, 0, 16000, 400, 0, 16000);
+                }
+                if (!hist2d1)
+                {
+                    hist2d1 = new TH2F(histName1, histName1,
+                                       800, -1, 1, 800, 0, 4000);
+                }
+
+                if (sx3EBk > 100 || sx3EUp > 100 || sx3EDn > 100)
+                {
+                    hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
+
+                    // Use the correlated triplet directly
+                    dataPoints[{detID, sx3ChBk, sx3ChUp, sx3ChDn}]
+                        .emplace_back(sx3EBk, sx3EUp, sx3EDn);
+                }
+
+                hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
+                hist2d1->Fill((sx3EUp - sx3EDn) / (sx3EUp + sx3EDn), sx3EBk);
+            }
+        }
+    }
+
+    return kTRUE;
+}
+
+const double GAIN_ACCEPTANCE_THRESHOLD = 0.3;
+void GainMatchSX3::Terminate()
+{
+    double backSlope[MAX_DET][MAX_BK] = {{0}};
+    bool backSlopeValid[MAX_DET][MAX_BK] = {{false}};
+
+    std::ofstream outFile("sx3_BackGains0.txt");
+    if (!outFile.is_open())
+    {
+        std::cerr << "Error opening sx3_BackGains.txt for writing!" << std::endl;
+        return;
+    }
+
+    // === Gain fit: (Up+Dn) vs Back, grouped by [id][bk] ===
+    for (int id = 0; id < MAX_DET; id++)
+    {
+        for (int bk = 0; bk < MAX_BK; bk++)
+        {
+            std::vector<double> bkE, udE;
+
+            // Collect all (Up+Dn, Back) for this id,bk
+            for (const auto &kv : dataPoints)
+            {
+                auto [cid, cbk, u, d] = kv.first;
+                if (cid != id || cbk != bk)
+                    continue;
+
+                for (const auto &pr : kv.second)
+                {
+                    double eBk, eUp, eDn;
+                    std::tie(eBk, eUp, eDn) = pr;
+                    if ((eBk < 100) || (eUp < 100) || (eDn < 100))
+                        continue;
+
+                    bkE.push_back(eBk);
+                    udE.push_back(eUp + eDn);
+                }
+            }
+
+            if (bkE.size() < 5)
+                continue; // not enough statistics
+
+            // Build graph with errors
+            const double fixedError = 0.0;                      // ADC channels
+            std::vector<double> exVals(udE.size(), 0.0);        // no x error
+            std::vector<double> eyVals(udE.size(), fixedError); // constant y error
+
+            TGraphErrors g(udE.size(), udE.data(), bkE.data(),
+                           exVals.data(), eyVals.data());
+
+            TF1 f("f", "[0]*x", 0, 16000);
+            // f.SetParameter(0, 1.0); // initial slope
+
+            if (drawCanvases)
+            {
+                g.SetTitle(Form("Detector %d Back %d: (Up+Dn) vs Back", id, bk));
+                g.SetMarkerStyle(20);
+                g.SetMarkerColor(kBlue);
+                g.Draw("AP");
+
+                g.Fit(&f, interactiveMode ? "Q" : "QNR");
+
+                if (verboseFit)
+                {
+                    double chi2 = f.GetChisquare();
+                    int ndf = f.GetNDF();
+                    double reducedChi2 = (ndf != 0) ? chi2 / ndf : -1;
+
+                    std::cout << Form("Det%d Back%d → Slope: %.4f | χ²/ndf = %.2f/%d = %.2f",
+                                      id, bk, f.GetParameter(0), chi2, ndf, reducedChi2)
+                              << std::endl;
+                }
+
+                if (interactiveMode)
+                {
+                    c.Update();
+                    gPad->WaitPrimitive();
+                }
+                else
+                {
+                    c.Close();
+                }
+            }
+            else
+            {
+                g.Fit(&f, "QNR");
+            }
+
+            double slope = 1 / f.GetParameter(0);
+            if (std::abs(slope - 1.0) < 0.3) // sanity check
+            {
+                backSlope[id][bk] = slope;
+                backSlopeValid[id][bk] = true;
+                outFile << id << " " << bk << " " << slope << "\n";
+                printf("Back slope Det%d Bk%d → %.4f\n", id, bk, slope);
+            }
+            else
+            {
+                std::cerr << "Warning: Bad slope for Det" << id << " Bk" << bk
+                          << " slope=" << slope << std::endl;
+            }
+        }
+    }
+
+    outFile.close();
+    std::cout << "Back gain matching complete." << std::endl;
+
+    // === Create histograms ===
+    TH2F *hFVB = new TH2F("hFVB", "Corrected Up+Dn vs Corrected Back;Up+Dn E;Corrected Back E",
+                          600, 0, 16000, 600, 0, 16000);
+    TH2F *hAsym = new TH2F("hAsym", "Up vs Dn divide corrected back;Up/Back E;Dn/Back E",
+                           400, 0.0, 1.0, 400, 0.0, 1.0);
+    TH2F *hAsymUnorm = new TH2F("hAsymUnorm", "Up vs Dn;Up E;Dn E",
+                                800, 0.0, 4000.0, 800, 0.0, 4000.0);
+
+    // Fill histograms using corrected back energies
+    for (const auto &kv : dataPoints)
+    {
+        auto [id, bk, u, d] = kv.first;
+        if (!backSlopeValid[id][bk])
+            continue;
+
+        double slope = backSlope[id][bk];
+
+        for (const auto &pr : kv.second)
+        {
+            double eBk, eUp, eDn;
+            std::tie(eBk, eUp, eDn) = pr;
+
+            double updn = eUp + eDn;
+            if (updn == 0 || eBk == 0)
+                continue;
+
+            double correctedBack = eBk * slope;
+            double asym = (eUp - eDn) / updn;
+
+            hFVB->Fill(updn, correctedBack);
+            hAsym->Fill(eUp / correctedBack, eDn / correctedBack);
+            hAsymUnorm->Fill(eUp, eDn);
+            TString histNamex = Form("CorrBack_id%d_U%d-D%dvsB%d", id, u, d, bk);
+
+            TH2F *hist2dx = (TH2F *)gDirectory->Get(histNamex);
+            if (!hist2dx)
+            {
+                hist2dx = new TH2F(histNamex, histNamex,
+                                   800, -1, 1, 800, 0, 4000);
+            }
+
+            hist2dx->Fill((eUp - eDn) / (eUp + eDn), correctedBack);
+        }
+    }
+    // plotter.FlushToDisk();
+}

GainMatchSX3.h

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

GainMatchSX3Front.C

@@ -0,0 +1,420 @@
+#define GainMatchSX3Front_cxx
+
+#include "GainMatchSX3Front.h"
+#include <TH2.h>
+#include <TF1.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TCutG.h>
+#include <fstream>
+#include <utility>
+#include <algorithm>
+#include <TProfile.h>
+#include "Armory/ClassSX3.h"
+#include "TGraphErrors.h"
+#include "TMultiDimFit.h"
+
+#include "TVector3.h"
+
+TH2F *hSX3FvsB;
+TH2F *hSX3FvsB_g;
+TH2F *hsx3IndexVE;
+TH2F *hsx3IndexVE_g;
+TH2F *hSX3;
+TH2F *hsx3Coin;
+
+int padID = 0;
+
+SX3 sx3_contr;
+TCutG *cut;
+TCutG *cut1;
+std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
+TCanvas c(Form("canvas"), "Fit", 800, 600);
+
+// Gain arrays
+
+const int MAX_DET = 24;
+const int MAX_UP = 4;
+const int MAX_DOWN = 4;
+const int MAX_BK = 4;
+double backGain[MAX_DET][MAX_BK] = {{0}};
+bool backGainValid[MAX_DET][MAX_BK] = {{false}};
+double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
+bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
+double uvdslope[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
+// ==== Configuration Flags ====
+const bool interactiveMode = true; // If true: show canvas + wait for user
+const bool verboseFit = true;      // If true: print fit summary and chi²
+const bool drawCanvases = true;    // If false: canvases won't be drawn at all
+
+void GainMatchSX3Front::Begin(TTree * /*tree*/)
+{
+    TString option = GetOption();
+
+    hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 800, 0, 16000, 800, 0, 16000);
+    hSX3FvsB_g = new TH2F("hSX3FvsB_g", "SX3 Front vs Back; Front E; Back E", 800, 0, 16000, 800, 0, 16000);
+    hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
+    hsx3IndexVE_g = new TH2F("hsx3IndexVE_g", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
+    hSX3 = new TH2F("hSX3", "SX3 Front v Back; Fronts; Backs", 8, 0, 8, 4, 0, 4);
+
+    hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12);
+
+    sx3_contr.ConstructGeo();
+
+    // Load the TCutG object
+    TFile *cutFile = TFile::Open("sx3cut.root");
+    bool cutLoaded = (cut != nullptr);
+    cut = dynamic_cast<TCutG *>(cutFile->Get("sx3cut"));
+    if (!cut)
+    {
+        std::cerr << "Error: Could not find TCutG named 'sx3cut' in sx3cut.root" << std::endl;
+        return;
+    }
+    cut->SetName("sx3cut"); // Ensure the cut has the correct name
+
+    // Load the TCutG object
+    TFile *cutFile1 = TFile::Open("UvD.root");
+    bool cut1Loaded = (cut1 != nullptr);
+    cut1 = dynamic_cast<TCutG *>(cutFile1->Get("UvD"));
+    if (!cut1)
+    {
+        std::cerr << "Error: Could not find TCutG named 'UvD' in UvD.root" << std::endl;
+        return;
+    }
+    cut1->SetName("UvD");
+
+    std::string filename = "sx3_BackGains.txt";
+
+    std::ifstream infile(filename);
+    if (!infile.is_open())
+    {
+        std::cerr << "Error opening " << filename << "!" << std::endl;
+        return;
+    }
+
+    int id, bk;
+    double gain;
+    while (infile >> id >> bk >> gain)
+    {
+        backGain[id][bk] = gain;
+        if (backGain[id][bk] > 0)
+            backGainValid[id][bk] = true;
+        else
+            backGainValid[id][bk] = false;
+    }
+
+    SX3 sx3_contr;
+}
+
+Bool_t GainMatchSX3Front::Process(Long64_t entry)
+{
+
+    b_sx3Multi->GetEntry(entry);
+    b_sx3ID->GetEntry(entry);
+    b_sx3Ch->GetEntry(entry);
+    b_sx3E->GetEntry(entry);
+    b_sx3T->GetEntry(entry);
+
+    sx3.CalIndex();
+
+    std::vector<std::pair<int, int>> ID;
+    for (int i = 0; i < sx3.multi; i++)
+    {
+
+        for (int j = i + 1; j < sx3.multi; j++)
+        {
+            // if (sx3.id[i] == 3)
+            hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
+        }
+        if (sx3.e[i] > 100)
+        {
+            ID.push_back(std::pair<int, int>(sx3.id[i], i));
+            hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
+        }
+    }
+
+    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; });
+
+        // start with the first entry in the sorted array: channels that belong to the same detector are together in sequenmce
+        std::vector<std::pair<int, int>> sx3ID;
+        sx3ID.push_back(ID[0]);
+        bool found = false;
+
+        for (size_t i = 1; i < ID.size(); i++)
+        { // Check if id of i belongs to the same detector and then add it to the detector ID vector
+            if (ID[i].first == sx3ID.back().first)
+            { // count the nunmber of hits that belong to the same detector
+                sx3ID.push_back(ID[i]);
+
+                if (sx3ID.size() >= 3)
+                {
+                    found = true;
+                }
+            }
+            else
+            { // the next event does not belong to the same detector, abandon the first event and continue with the next one
+                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;
+                // Check the channel number and assign it to the appropriate channel type
+                if (sx3.ch[index] < 8)
+                {
+                    if (sx3.ch[index] % 2 == 0)
+                    {
+                        sx3ChDn = sx3.ch[index] / 2;
+                        sx3EDn = sx3.e[index];
+                    }
+                    else
+                    {
+                        sx3ChUp = sx3.ch[index] / 2;
+                        sx3EUp = sx3.e[index];
+                    }
+                }
+                else
+                {
+                    sx3ChBk = sx3.ch[index] - 8;
+                    // if (sx3ChBk == 2)
+                    //     printf("Found back channel Det %d Back %d \n", sx3.id[index], sx3ChBk);
+                    sx3EBk = sx3.e[index];
+                }
+            }
+
+            for (int i = 0; i < sx3.multi; i++)
+            {
+                // If we have a valid front and back channel, fill the histograms
+                hSX3->Fill(sx3ChDn + 4, sx3ChBk);
+                hSX3->Fill(sx3ChUp, sx3ChBk);
+
+                // Fill the histogram for the front vs back
+                hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
+
+                if (sx3.e[i] > 100 && sx3.id[i] == 3)
+                {
+                    // back gain correction
+
+                    // Fill the histogram for the front vs back with gain correction
+                    // hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
+                    // // Fill the index vs energy histogram
+                    // hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
+                    // }
+                    // {
+                    TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk);
+                    TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
+                    if (!hist2d)
+                    {
+                        hist2d = new TH2F(histName, Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk), 400, 0, 16000, 400, 0, 16000);
+                    }
+
+                    // if (sx3ChBk == 2)
+                    //     printf("Found back channel Det %d Back %d \n", sx3.id[i], sx3ChBk);
+                    hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
+                    hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
+
+                    hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
+
+                    if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk) && cut1 && cut1->IsInside(sx3EUp / sx3EBk, sx3EDn / sx3EBk))
+                    {
+
+                        if (backGainValid[sx3.id[i]][sx3ChBk])
+                        {
+                            sx3EBk *= backGain[sx3.id[i]][sx3ChBk];
+                        }
+                        // Accumulate data for gain matching
+                        dataPoints[{sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(sx3EBk, sx3EUp, sx3EDn);
+                    }
+                }
+            }
+        }
+    }
+
+    return kTRUE;
+}
+
+void GainMatchSX3Front::Terminate()
+{
+
+    std::map<std::tuple<int, int, int, int>, TVectorD> fitCoefficients;
+
+    // === Gain matching ===
+
+    std::ofstream outFile("sx3_GainMatchfront.txt");
+    if (!outFile.is_open())
+    {
+        std::cerr << "Error opening output file!" << std::endl;
+        return;
+    }
+
+    TH2F *hUvD = new TH2F("hUvD", " UvD; Up/CorrBack; Down/CorrBack", 600, 0, 1, 600, 0, 1);
+
+    for (const auto &kv : dataPoints)
+    {
+        auto [id, bk, u, d] = kv.first;
+        const auto &pts = kv.second;
+
+        if (pts.size() < 50)
+            continue;
+
+        std::vector<double> uE, dE, udE, corrBkE;
+
+        for (const auto &pr : pts)
+        {
+            double eBkCorr, eUp, eDn;
+            std::tie(eBkCorr, eUp, eDn) = pr;
+            if ((eBkCorr < 100) || (eUp < 100) || (eDn < 100))
+                continue; // Skip if any energy is zero
+            uE.push_back(eUp / eBkCorr);
+            dE.push_back(eDn / eBkCorr);
+            udE.push_back(eUp + eDn);
+            corrBkE.push_back(eBkCorr);
+            hUvD->Fill(eUp / eBkCorr, eDn / eBkCorr);
+        }
+        if (uE.size() < 5 || dE.size() < 5 || corrBkE.size() < 5)
+            continue; // Ensure we have enough points for fitting
+        // TGraph g(udE.size(), udE.data(), corrBkE.data());
+
+        // TF1 f("f", "[0]*x", 0, 20000);
+        // f.SetParameter(0, 1.0); // Initial guess for the gain
+        // g.Fit(&f, "R");
+
+        const double fixedError = 0.0; // in ADC channels
+
+        std::vector<double> xVals, yVals, exVals, eyVals;
+
+        // Build data with fixed error
+        for (size_t i = 0; i < udE.size(); ++i)
+        {
+            double x = uE[i]; // front energy
+            double y = dE[i]; // back energy
+
+            xVals.push_back(x);
+            yVals.push_back(y);
+            exVals.push_back(fixedError); // error in up energy
+            eyVals.push_back(0.);         // error in down energy
+        }
+
+        // Build TGraphErrors with errors
+        TGraphErrors g(xVals.size(), xVals.data(), yVals.data(), exVals.data(), eyVals.data());
+
+        TF1 f("f", "[0]*x+[1]", 0, 16000);
+        f.SetParameter(0, -1.0); // Initial guess
+
+        if (drawCanvases)
+        {
+            g.SetTitle(Form("Detector %d: U%d D%d B%d", id, u, d, bk));
+            g.SetMarkerStyle(20);
+            g.SetMarkerColor(kBlue);
+            g.Draw("AP");
+
+            g.Fit(&f, interactiveMode ? "Q" : "QNR"); // 'R' avoids refit, 'N' skips drawing
+
+            if (verboseFit)
+            {
+                double chi2 = f.GetChisquare();
+                int ndf = f.GetNDF();
+                double reducedChi2 = (ndf != 0) ? chi2 / ndf : -1;
+
+                std::cout << Form("Det%d U%d D%d B%d → Gain: %.4f | χ²/ndf = %.2f/%d = %.2f",
+                                  id, u, d, bk, f.GetParameter(0), chi2, ndf, reducedChi2)
+                          << std::endl;
+            }
+
+            if (interactiveMode)
+            {
+                c.Update();
+                gPad->WaitPrimitive();
+            }
+            else
+            {
+                c.Close(); // Optionally avoid clutter in batch
+            }
+        }
+        else
+        {
+            g.Fit(&f, "QNR");
+        }
+
+        double slope = f.GetParameter(0);
+        double intercept = f.GetParameter(1);   
+
+        // printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]);
+        if (std::abs(slope + 1.0) < 0.3) // sanity check
+        {
+            frontGain[id][bk][u][d] = slope;
+
+            frontGainValid[id][bk][u][d] = true;
+            outFile << id << " " << bk << " " << u << " " << d << " " << TMath::Abs(slope)/intercept << " " << 1.0/intercept << std::endl;
+            printf("Back slope Det%d Bk%d → %.4f\n", id, bk, slope);
+        }
+        else
+        {
+            std::cerr << "Warning: Bad slope for Det" << id << " Bk" << bk
+                      << " slope=" << f.GetParameter(0) << std::endl;
+        }
+    }
+
+    outFile.close();
+    std::cout << "Gain matching complete." << std::endl;
+
+    // === Stage 3: Create corrected histogram ===
+    TH2F *hCorrectedFvB = new TH2F("hCorrectedFvB", "Corrected;Corrected Front Sum;Corrected Back", 800, 0, 8000, 800, 0, 8000);
+    TH2F *hCorrectedUvD = new TH2F("hCorrectedUvD", "Corrected UvD; UvD Up; UvD Down", 600, 0, 1, 600, 0, 1);
+
+    for (const auto &kv : dataPoints)
+    {
+
+        auto [id, bk, u, d] = kv.first;
+        double front;
+        if (frontGainValid[id][bk][u][d])
+            front = frontGain[id][bk][u][d];
+        else
+            continue;
+        for (const auto &pr : kv.second)
+        {
+            double eBk, eUp, eDn;
+            std::tie(eBk, eUp, eDn) = pr;
+            double corrUp = eUp * front;
+            // double corrDn = eDn * front;
+
+            hCorrectedFvB->Fill(corrUp + eDn, eBk);
+            hCorrectedUvD->Fill(corrUp / eBk, eDn / eBk);
+        }
+    }
+
+    // // === Final canvas ===
+    // gStyle->SetOptStat(1110);
+    // TCanvas *c1 = new TCanvas("c1", "Gain Correction Results", 1200, 600);
+    // c1->Divide(2, 1);
+
+    // c1->cd(1);
+    // hSX3FvsB_g->SetTitle("Before Correction (Gated)");
+    // hSX3FvsB_g->GetXaxis()->SetTitle("Measured Front Sum (E_Up + E_Dn)");
+    // hSX3FvsB_g->GetYaxis()->SetTitle("Measured Back E");
+    // hSX3FvsB_g->Draw("colz");
+
+    // c1->cd(2);
+    // hCorrectedFvB->SetTitle("After Correction");
+    // hCorrectedFvB->Draw("colz");
+    // TF1 *diag = new TF1("diag", "x", 0, 40000);
+    // diag->SetLineColor(kRed);
+    // diag->SetLineWidth(2);
+    // diag->Draw("same");
+
+    std::cout << "Terminate() completed successfully." << std::endl;
+}

GainMatchSX3Front.h

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

GainMatchSX3Front1.C

@@ -0,0 +1,245 @@
+#define GainMatchSX3_cxx
+
+#include "GainMatchSX3.h"
+#include "Armory/ClassSX3.h"
+#include <TFile.h>
+#include <TTree.h>
+#include <TGraph.h>
+#include <TF1.h>
+#include <TH2F.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TApplication.h>
+#include <map>
+#include <vector>
+#include <tuple>
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+// Constants
+const int MAX_DET = 24;
+const int MAX_BK = 4;
+const int MAX_UP = 4;
+const int MAX_DOWN = 4;
+
+// Gain arrays
+double backGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
+bool backGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
+
+double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
+bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
+
+// Data container
+std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
+
+// Load back gains
+void LoadBackGains(const std::string &filename)
+{
+    std::ifstream infile(filename);
+    if (!infile.is_open())
+    {
+        std::cerr << "Error opening " << filename << "!" << std::endl;
+        return;
+    }
+
+    int id, bk, u, d;
+    double gain;
+    while (infile >> id >> bk >> u >> d >> gain)
+    {
+        backGain[id][bk][u][d] = gain;
+        backGainValid[id][bk][u][d] = true;
+    }
+
+    infile.close();
+    std::cout << "Loaded back gains from " << filename << std::endl;
+    SX3 sx3_contr;
+}
+
+// Front gain matching function
+Bool_t GainMatchSX3::Process(Long64_t entry)
+{
+    // Link SX3 branches
+    
+    b_sx3Multi->GetEntry(entry);
+    b_sx3ID->GetEntry(entry);
+    b_sx3Ch->GetEntry(entry);
+    b_sx3E->GetEntry(entry);
+    b_sx3T->GetEntry(entry);
+
+    sx3.CalIndex();
+
+    Long64_t nentries = tree->GetEntries(Long64_t entry);
+    std::cout << "Total entries: " << nentries << std::endl;
+
+    TH2F *hBefore = new TH2F("hBefore", "Before Correction;E_Up+E_Dn;Back Energy", 400, 0, 40000, 400, 0, 40000);
+    TH2F *hAfter = new TH2F("hAfter", "After Correction;E_Up+E_Dn;Corrected Back Energy", 400, 0, 40000, 400, 0, 40000);
+
+    for (Long64_t entry = 0; entry < nentries; ++entry)
+    {
+        tree->GetEntry(entry);
+        sx3.CalIndex();
+
+        std::vector<std::pair<int, int>> ID;
+
+        for (int i = 0; i < sx3.multi; i++)
+        {
+            if (sx3.e[i] > 100)
+            {
+                ID.push_back({sx3.id[i], i});
+            }
+        }
+
+        if (ID.empty())
+            continue;
+
+        // Sort by id
+        std::sort(ID.begin(), ID.end(), [](auto &a, auto &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)
+            continue;
+
+        int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1;
+        float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0;
+        int sx3id = sx3ID[0].first;
+
+        for (auto &[id, idx] : sx3ID)
+        {
+            if (sx3.ch[idx] < 8)
+            {
+                if (sx3.ch[idx] % 2 == 0)
+                {
+                    sx3ChDn = sx3.ch[idx] / 2;
+                    sx3EDn = sx3.e[idx];
+                }
+                else
+                {
+                    sx3ChUp = sx3.ch[idx] / 2;
+                    sx3EUp = sx3.e[idx];
+                }
+            }
+            else
+            {
+                sx3ChBk = sx3.ch[idx] - 8;
+                sx3EBk = sx3.e[idx];
+            }
+        }
+
+        if (sx3ChUp < 0 || sx3ChDn < 0 || sx3ChBk < 0)
+            continue;
+
+        if (!backGainValid[sx3id][sx3ChBk][sx3ChUp][sx3ChDn])
+            continue;
+
+        double corrBk = sx3EBk * backGain[sx3id][sx3ChBk][sx3ChUp][sx3ChDn];
+
+        hBefore->Fill(sx3EUp + sx3EDn, sx3EBk);
+        hAfter->Fill(sx3EUp + sx3EDn, corrBk);
+
+        dataPoints[{sx3id, sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(corrBk, sx3EUp, sx3EDn);
+    }
+
+    // === Fit front gains ===
+    std::ofstream outFile("sx3_GainMatchfront.txt");
+    if (!outFile.is_open())
+    {
+        std::cerr << "Error opening sx3_GainMatchfront.txt!" << std::endl;
+        return;
+    }
+
+    for (const auto &kv : dataPoints)
+    {
+        auto [id, bk, u, d] = kv.first;
+        const auto &pts = kv.second;
+
+        if (pts.size() < 5)
+            continue;
+
+        std::vector<double> udE, corrBkE;
+
+        for (const auto &pr : pts)
+        {
+            double eBkCorr, eUp, eDn;
+            std::tie(eBkCorr, eUp, eDn) = pr;
+            udE.push_back(eUp + eDn);
+            corrBkE.push_back(eBkCorr);
+        }
+
+        TGraph g(udE.size(), udE.data(), corrBkE.data());
+        TF1 f("f", "[0]*x", 0, 40000);
+        g.Fit(&f, "QNR");
+
+        frontGain[id][bk][u][d] = f.GetParameter(0);
+        frontGainValid[id][bk][u][d] = true;
+
+        outFile << id << " " << bk << " " << u << " " << d << " " << frontGain[id][bk][u][d] << std::endl;
+        printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]);
+    }
+
+    outFile.close();
+    std::cout << "Front gain matching complete." << std::endl;
+
+    // === Draw diagnostic plots ===
+    gStyle->SetOptStat(1110);
+    TCanvas *c = new TCanvas("c", "Gain Matching Diagnostics", 1200, 600);
+    c->Divide(2, 1);
+
+    c->cd(1);
+    hBefore->Draw("colz");
+    TF1 *diag1 = new TF1("diag1", "x", 0, 40000);
+    diag1->SetLineColor(kRed);
+    diag1->Draw("same");
+
+    c->cd(2);
+    hAfter->Draw("colz");
+    TF1 *diag2 = new TF1("diag2", "x", 0, 40000);
+    diag2->SetLineColor(kRed);
+    diag2->Draw("same");
+}
+
+int main(int argc, char **argv)
+{
+    TApplication app("app", &argc, argv);
+
+    // Load back gains
+    LoadBackGains("sx3_GainMatchback.txt");
+
+    // Open tree
+    TFile *f = TFile::Open("input_tree.root"); // <<< Change file name
+    if (!f || f->IsZombie())
+    {
+        std::cerr << "Cannot open input_tree.root!" << std::endl;
+        return 1;
+    }
+    TTree *tree = (TTree *)f->Get("tree");
+    if (!tree)
+    {
+        std::cerr << "Tree not found!" << std::endl;
+        return 1;
+    }
+
+    // Run front gain matching
+    GainMatchSX3(tree);
+
+    app.Run();
+    return 0;
+}

MatchAndPlotCentroids.C

@@ -0,0 +1,133 @@
+#include <fstream>
+#include <sstream>
+#include <vector>
+#include <map>
+#include <iostream>
+#include <TGraph.h>
+#include <TF1.h>
+#include <TCanvas.h>
+#include <TH1.h>
+
+
+
+void MatchAndPlotCentroids() {
+    // Open the centroid data file
+    std::ifstream inputFile("centroids.txt");
+    if (!inputFile.is_open()) {
+        std::cerr << "Error: Could not open Centroids.txt" << std::endl;
+        return;
+    }
+
+    // Data structure to store centroids by histogram and peak number
+    std::map<int, std::map<int, double>> centroidData;
+
+    // Read data from the file
+    std::string line;
+    while (std::getline(inputFile, line)) {
+        std::istringstream iss(line);
+        int histogramIndex, peakNumber;
+        double centroid;
+        if (iss >> histogramIndex >> peakNumber >> centroid) {
+            centroidData[histogramIndex][peakNumber] = centroid;
+        }
+    }
+
+    inputFile.close();
+
+    // Ensure histogram 24 exists and has data
+    if (centroidData.find(1) == centroidData.end()) {
+        std::cerr << "Error: Histogram 0 not found in the data!" << std::endl;
+        return;
+    }
+
+    // Reference centroids from histogram 24
+    const auto& referenceCentroids = centroidData[1];
+    std::ofstream outputFile("slope_intercept_results.txt");
+    if (!outputFile.is_open()) {
+        std::cerr << "Error: Could not open the output file for writing!" << std::endl;
+        return;
+    }
+    outputFile << "Histogram Number\tSlope\tIntercept\n";
+    // Loop through histograms 25 to 47
+    for (int targetHist = 0; targetHist <= 23; targetHist++) {
+        // Ensure the target histogram exists and matches in peak numbers
+        if (centroidData.find(targetHist) == centroidData.end() || centroidData[targetHist].size() != referenceCentroids.size()) {
+            //4th cnetroid data point for 19 was generated using the 3 datqa points for the slope of wires 0 and 19
+            std::cout << "Skipping Histogram " << targetHist << " due to mismatched or missing data." << std::endl;
+            continue;
+        }
+
+        // Prepare x and y values for TGraph
+        std::vector<double> xValues, yValues;
+        for (const auto& [peakNumber, refCentroid] : referenceCentroids) {
+            if (centroidData[targetHist].find(peakNumber) != centroidData[targetHist].end()) {
+                yValues.push_back(refCentroid);
+                xValues.push_back(centroidData[targetHist][peakNumber]);
+            } else {
+                std::cerr << "Warning: Peak " << peakNumber << " missing in histogram " << targetHist << std::endl;
+            }
+        }
+
+        if (xValues.size() < 3) {
+            std::cout << "Skipping Histogram " << targetHist << " as it has less than 3 matching centroids." << std::endl;
+            continue;
+        }
+
+        // Create a TGraph
+        TCanvas *c1 = new TCanvas(Form("c_centroid_1_vs_%d", targetHist), Form("Centroid 1 vs %d", targetHist), 800, 600);
+        TGraph *graph = new TGraph(xValues.size(), &xValues[0], &yValues[0]);
+        graph->SetTitle(Form("Centroid of Histogram  %d vs 1", targetHist));
+        graph->GetYaxis()->SetTitle("Centroid of Histogram 1");
+        graph->GetXaxis()->SetTitle(Form("Centroid of Histogram %d", targetHist));
+        graph->SetMarkerStyle(20); // Full circle marker
+        graph->SetMarkerSize(1.0);
+        graph->SetMarkerColor(kBlue);
+        // Draw the graph
+        graph->Draw("AP");
+        double minX = *std::min_element(xValues.begin(), xValues.end());
+        double maxX = *std::max_element(xValues.begin(), xValues.end());
+        // Fit the data with a linear function
+        TF1 *fitLine = new TF1("fitLine", "pol1", minX, maxX);  // Adjust range as needed
+	fitLine->SetLineColor(kRed); // Set the line color to distinguish it
+        fitLine->SetLineWidth(2);    // Thicker line for visibility
+        graph->Fit(fitLine, "M");
+        fitLine->Draw("same");
+        fitLine->SetParLimits(0, -10, 10); // Limit intercept between -10 and 10
+        fitLine->SetParLimits(1, 0, 2); 
+        // Extract slope and intercept
+        double slope = fitLine->GetParameter(1);
+        double intercept = fitLine->GetParameter(0);
+        outputFile << targetHist << "\t" << slope << "\t" << intercept << "\n";
+        std::cout << "Histogram 24 vs " << targetHist << ": Slope = " << slope << ", Intercept = " << intercept << std::endl;
+        std::vector<double> residuals;
+    for (size_t i = 0; i < xValues.size(); ++i) {
+        double fittedY = fitLine->Eval(xValues[i]); // Evaluate fitted function at x
+        double residual = yValues[i] - fittedY;    // Residual = observed - fitted
+        residuals.push_back(residual);
+    }
+
+    // Create a graph for the residuals
+    /*TGraph *residualGraph = new TGraph(residuals.size(), &xValues[0], &residuals[0]);
+    residualGraph->SetTitle(Form("Residuals for Histogram 24 vs %d", targetHist));
+    residualGraph->GetYaxis()->SetTitle("Residuals");
+    residualGraph->GetXaxis()->SetTitle(Form("Centroid of Histogram %d", targetHist));
+    residualGraph->SetMarkerStyle(20);
+    residualGraph->SetMarkerSize(1.0);
+    residualGraph->SetMarkerColor(kGreen);
+
+    // Draw the residuals plot below the original plot (can be on a new canvas if preferred)
+    TCanvas *c2 = new TCanvas(Form("c_residuals_24_vs_%d", targetHist), Form("Residuals for Centroid 24 vs %d", targetHist), 800, 400);
+    residualGraph->Draw("AP");*/
+        c1->Update();
+        //c2->Update();
+        std::cout << "Press Enter to continue..." << std::endl;
+        
+        //std::cin.get();
+        c1->WaitPrimitive();
+        //c2->WaitPrimitive();
+        //std::cin.get();
+        //std::cin.get();
+    }
+    outputFile.close();
+    std::cout << "Results written to slope_intercept_results.txt" << std::endl;
+}

PCGainMatch.C

@@ -0,0 +1,454 @@
+#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();
+}

PCGainMatch.h

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

ProcessRun.sh

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

QQQ_Calcheck.C

@@ -0,0 +1,167 @@
+
+#define QQQ_Calcheck_cxx
+
+#include <TH2.h>
+#include <TF1.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TCutG.h>
+#include <fstream>
+#include <utility>
+#include <algorithm>
+#include "Armory/HistPlotter.h"
+#include "TVector3.h"
+#include "QQQ_Calcheck.h"
+
+TH2F *hQQQFVB;
+HistPlotter *plotter;
+int padID = 0;
+
+TCutG *cut;
+std::map<std::tuple<int, int, int>, std::vector<std::pair<double, double>>> dataPoints;
+
+bool qqqEcut = false;
+
+// Gain Arrays
+const int MAX_QQQ = 4;
+const int MAX_RING = 16;
+const int MAX_WEDGE = 16;
+double qqqwGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
+double qqqrGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
+bool qqqwGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
+bool qqqrGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
+double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
+bool qqqCalibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
+
+void QQQ_Calcheck::Begin(TTree * /*tree*/)
+{
+    plotter = new HistPlotter("Cal_checkQQQ.root", "TFILE");
+    // ----------------------- Load QQQ Gains
+    {
+        std::string filename = "qqq_GainMatch.txt";
+        std::ifstream infile(filename);
+        if (!infile.is_open())
+        {
+            std::cerr << "Error opening " << filename << "!" << std::endl;
+        }
+        else
+        {
+            int det, ring, wedge;
+            double gainw,gainr;
+            while (infile >> det >> ring >> wedge >> gainw>> gainr)
+            {
+                qqqwGain[det][ring][wedge] = gainw;
+                // qqqrGain[det][ring][wedge] = gainr;
+                qqqwGainValid[det][ring][wedge] = (gainw > 0);
+                // qqqrGainValid[det][ring][wedge] = (gainr > 0);
+            }
+            infile.close();
+            std::cout << "Loaded QQQ gains from " << filename << std::endl;
+        }
+    }
+    // ----------------------- Load QQQ Calibrations
+    {
+        std::string filename = "qqq_Calib.txt";
+        std::ifstream infile(filename);
+        if (!infile.is_open())
+        {
+            std::cerr << "Error opening " << filename << "!" << std::endl;
+        }
+        else
+        {
+            int det, ring, wedge;
+            double slope;
+            while (infile >> det >> ring >> wedge >> slope)
+            {
+                qqqCalib[det][ring][wedge] = slope;
+                qqqCalibValid[det][ring][wedge] = (slope > 0);
+            }
+            infile.close();
+            std::cout << "Loaded QQQ calibrations from " << filename << std::endl;
+        }
+    }
+}
+
+Bool_t QQQ_Calcheck::Process(Long64_t entry)
+{
+    b_qqqMulti->GetEntry(entry);
+    b_qqqID->GetEntry(entry);
+    b_qqqCh->GetEntry(entry);
+    b_qqqE->GetEntry(entry);
+    b_qqqT->GetEntry(entry);
+
+    qqq.CalIndex();
+
+    for (int i = 0; i < qqq.multi; i++)
+    {
+        for (int j = i + 1; j < qqq.multi; j++)
+        {
+            if (qqq.e[i] > 100)
+                qqqEcut = true;
+            if (qqq.id[i] == qqq.id[j])
+            {
+                int chWedge = -1;
+                int chRing = -1;
+                float eWedgeRaw = 0.0;
+                float eWedge = 0.0;
+                float eWedgeMeV = 0.0;
+                float eRingRaw = 0.0;
+                float eRing = 0.0;
+                float eRingMeV = 0.0;
+                // plug in gains
+                if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && /*qqqrGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16] &&*/ qqqwGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
+                {
+                    chWedge = qqq.ch[i];
+                    eWedgeRaw = qqq.e[i];
+                    eWedge = qqq.e[i] * qqqwGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
+                    // printf("Wedge E: %.2f  Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]);
+                    chRing = qqq.ch[j] - 16;
+                    eRingRaw = qqq.e[j];
+                    eRing = qqq.e[j];//* qqqrGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i]-16];
+                }
+                else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16/* && qqqrGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16] */&& qqqwGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
+                {
+                    chWedge = qqq.ch[j];
+                    eWedge = qqq.e[j] * qqqwGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
+                    eWedgeRaw = qqq.e[j];
+
+                    chRing = qqq.ch[i] - 16;
+                    eRing = qqq.e[i];// * qqqrGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
+                    eRingRaw = qqq.e[i];
+                }
+                else
+                    continue;
+                // plug in calibrations
+                if (qqqCalibValid[qqq.id[i]][chRing][chWedge])
+                {
+                    eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
+                    eRingMeV = eRing * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
+                }
+                else
+                    continue;
+
+                // hQQQFVB->Fill(eWedge, eRing);
+                plotter->Fill2D(Form("hRaw_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 8000, 400, 0, 8000, eWedgeRaw, eRingRaw, "ERaw");
+                plotter->Fill2D(Form("hGM_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000, eWedge, eRing, "EGM");
+                plotter->Fill2D(Form("hCal_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 10, 400, 0, 10, eWedgeMeV, eRingMeV, "ECal");
+                if(eWedgeRaw >1500 && eRingRaw>1500 )
+                plotter->Fill2D(Form("hCal_cut_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 10, 400, 0, 10, eWedgeMeV, eRingMeV, "ECal_cut");
+                
+                plotter->Fill2D(Form("hRCal_qqq%d", qqq.id[i]), 16, 0, 15, 1000, 0, 30, chRing, eRingMeV, "RingCal");
+                plotter->Fill2D(Form("hWCal_qqq%d", qqq.id[i]), 16, 0, 15, 1000, 0, 30, chWedge, eWedgeMeV, "WedgeCal");
+                plotter->Fill2D("hRawQQQ", 4000, 0, 8000, 4000, 0, 8000, eWedgeRaw, eRingRaw);
+                plotter->Fill2D("hGMQQQ", 4000, 0, 8000, 4000, 0, 8000, eWedge, eRing);
+                plotter->Fill2D("hCalQQQ", 4000, 0, 10, 4000, 0, 10, eWedgeMeV, eRingMeV);
+            }
+        }
+    }
+
+    return kTRUE;
+}
+
+void QQQ_Calcheck::Terminate()
+{
+    plotter->FlushToDisk();
+    std::cout << "Calibration check file for 2D QQQ histogram saved.\n";
+}

QQQ_Calcheck.h

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

centroids.txt

@@ -0,0 +1,97 @@
+HistogramIndex	PeakNumber	Centroid	Amplitude	Sigma
+0	1	991.118
+0	2	2026.83
+0	3	3060.26
+0	4	4092.45
+1	1	922.213
+1	2	1885.55
+1	3	2845.53
+1	4	3810.32
+2	1	955.591
+2	2	1953.17
+2	3	2949.37
+2	4	3950.79
+3	1	995.787
+3	2	2036.58
+3	3	3076.91
+3	4	4112.05
+4	1	1017.48
+4	2	2080.19
+4	3	3142.24
+4	4	4206.1
+5	1	1022.78
+5	2	2091.21
+5	3	3158.28
+5	4	4226.97
+6	1	1076.22
+6	2	2203.37
+6	3	3329.53
+6	4	4457.69
+7	1	977.46
+7	2	1998.02
+7	3	3017.36
+7	4	4040.47
+8	1	1049.74
+8	2	2144.38
+8	3	3238.2
+8	4	4335.25
+9	1	1000.59
+9	2	2046.42
+9	3	3090.29
+9	4	4129.63
+10	1	1014.92
+10	2	2076.16
+10	3	3134.59
+10	4	4213.42
+11	1	1004.85
+11	2	2052.88
+11	3	3100.3
+11	4	4164.75
+12	1	945.861
+12	2	1932.49
+12	3	2917.95
+12	4	3955.15
+13	1	998.307
+13	2	2040.38
+13	3	3078.76
+13	4	4135.51
+14	1	966.429
+14	2	1972.15
+14	3	2974.84
+14	4	4056.41
+15	1	958.352
+15	2	1958.64
+15	3	2957.7
+15	4	3970.41
+16	1	970.732
+16	2	1977.63
+16	3	2984.97
+16	4	4002.56
+17	1	1013.65
+17	2	2064.9
+17	3	3114.19
+17	4	4190.98
+18	1	975.538
+18	2	1990.64
+18	3	3005.46
+18	4	4048.99
+19	1	1082.91
+19	2	2194.08
+19	3	3303.65
+19	4	4411.32
+20	1	912.778
+20	2	1866.83
+20	3	2819.21
+20	4	3781.63
+21	1	1002.36
+21	2	1989.95
+21	3	2975.53
+21	4	3986.71
+22	1	1075.38
+22	2	2144.25
+22	3	3210.17
+22	4	4312.84
+23	1	988.828
+23	2	2016.35
+23	3	3044.19
+23	4	4082.41

centroids_edited.txt

@@ -0,0 +1,89 @@
+HistogramIndex	PeakNumber	Centroid	Amplitude	Sigma
+1	1	922.213
+1	2	1885.55
+1	3	2845.53
+1	4	3810.32
+2	1	955.591
+2	2	1953.17
+2	3	2949.37
+2	4	3950.79
+3	1	995.787
+3	2	2036.58
+3	3	3076.91
+3	4	4112.05
+4	1	1017.48
+4	2	2080.19
+4	3	3142.24
+4	4	4206.1
+5	1	1022.78
+5	2	2091.21
+5	3	3158.28
+5	4	4226.97
+6	1	1076.22
+6	2	2203.37
+6	3	3329.53
+6	4	4457.69
+7	1	977.46
+7	2	1998.02
+7	3	3017.36
+7	4	4040.47
+8	1	1049.74
+8	2	2144.38
+8	3	3238.2
+8	4	4335.25
+9	1	1000.59
+9	2	2046.42
+9	3	3090.29
+9	4	4129.63
+10	1	1014.92
+10	2	2076.16
+10	3	3134.59
+10	4	4213.42
+11	1	1004.85
+11	2	2052.88
+11	3	3100.3
+11	4	4164.75
+12	1	945.861
+12	2	1932.49
+12	3	2917.95
+12	4	3955.15
+13	1	998.307
+13	2	2040.38
+13	3	3078.76
+13	4	4135.51
+14	1	966.429
+14	2	1972.15
+14	3	2974.84
+14	4	4056.41
+15	1	958.352
+15	2	1958.64
+15	3	2957.7
+15	4	3970.41
+16	1	970.732
+16	2	1977.63
+16	3	2984.97
+16	4	4002.56
+17	1	1013.65
+17	2	2064.9
+17	3	3114.19
+17	4	4190.98
+18	1	975.538
+18	2	1990.64
+18	3	3005.46
+18	4	4048.99
+20	1	912.778
+20	2	1866.83
+20	3	2819.21
+20	4	3781.63
+21	1	1002.36
+21	2	1989.95
+21	3	2975.53
+21	4	3986.71
+22	1	1075.38
+22	2	2144.25
+22	3	3210.17
+22	4	4312.84
+23	1	988.828
+23	2	2016.35
+23	3	3044.19
+23	4	4082.41

makeplots.C

@@ -0,0 +1,87 @@
+#include "GainMatchQQQ.h"
+#include <TH2.h>
+#include <TF1.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TCutG.h>
+#include <fstream>
+#include <utility>
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+#include "Armory/HistPlotter.h"
+#include "TVector3.h"
+
+void makeplots()
+{
+    TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 1024, 768);
+    c1->Divide(2, 2);
+
+    TH2F *h1, *h2, *h3, *h4 = nullptr;
+
+    std::string type = "Cal_cut";
+    TFile *inFile = TFile::Open("Cal_checkQQQ.root");
+    TF1 f1("x","x",0 ,80);
+
+
+    for (int ring = 0; ring < 16; ring++)
+    {
+        for (int wedge = 0; wedge < 16; wedge++)
+        {
+
+            c1->cd(1);
+            h1 = (TH2F *)inFile->Get(Form("E%s/h%s_qqq%d_ring%d_wedge%d", type.c_str(), type.c_str(), 0, ring, wedge));
+            gPad->SetGrid();
+            if (h1)
+            {
+                h1->SetTitle(Form("QQQ%d Ring %d Wedge %d %s E", 0, ring, wedge, type.c_str()));
+                h1->Draw("COLZ");
+                f1.DrawClone("same"); 
+                
+                gPad->Modified();
+                gPad->Update();
+            }
+            
+            c1->cd(2);
+            h2 = (TH2F *)inFile->Get(Form("E%s/h%s_qqq%d_ring%d_wedge%d", type.c_str(), type.c_str(), 1, ring, wedge));
+            gPad->SetGrid();
+            if (h2)
+            {
+                std::cout<<h2<<std::endl;
+                h2->SetTitle(Form("QQQ%d Ring %d Wedge %d %s E", 1, ring, wedge, type.c_str()));
+                h2->Draw("COLZ");
+                f1.DrawClone("same");
+                gPad->Modified();
+                gPad->Update();
+            }
+            
+            c1->cd(3);
+            h3 = (TH2F *)inFile->Get(Form("E%s/h%s_qqq%d_ring%d_wedge%d", type.c_str(), type.c_str(), 2, ring, wedge));
+            gPad->SetGrid();
+            if (h3)
+            {
+                h3->SetTitle(Form("QQQ%d Ring %d Wedge %d %s E", 2, ring, wedge, type.c_str()));
+                h3->Draw("COLZ");
+                f1.DrawClone("same");
+                gPad->Modified();
+                gPad->Update();
+            }
+            
+            c1->cd(4);
+            h4 = (TH2F *)inFile->Get(Form("E%s/h%s_qqq%d_ring%d_wedge%d", type.c_str(), type.c_str(), 3, ring, wedge));
+            gPad->SetGrid();
+            if (h4)
+            {
+                h4->SetTitle(Form("QQQ%d Ring %d Wedge %d %s E", 3, ring, wedge, type.c_str()));
+                h4->Draw("COLZ");
+                f1.DrawClone("same");
+                gPad->Modified();
+                gPad->Update();
+            }
+            while (gPad->WaitPrimitive());
+        }
+    }
+
+    inFile->Close();
+}

mapping.h

@@ -19,6 +19,7 @@ const std::map<int, unsigned short> board = {
    {4, 22129},
    {5, 15529},
    {6, 15528},
+  //  {7,89},
    {7,   334},
    {8,   379},
    {9,   325},
@@ -27,13 +28,14 @@ const std::map<int, unsigned short> board = {
 const int nBd = board.size();
 
 const int nV1740 = 7;
-const int nV1725 = 3;
+const int nV1725 = 4;
 
 //+++++++++++++++++++ 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
 
@@ -76,17 +78,23 @@ 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, 20120, 20121, 20122, 20123, 20016, 20017, 20018, 20019, 20020, 20021, 20022, 20023,
+ 20116, 20117, 20118, 20119,    -1, 20121, 20122, 20123, 20016, 20017, 20018,    -1, 20020, 20021, 20022, 20023,
   //================== 379
- 20000, 20001, 20002, 20003, 20004, 20005,    -1, 20007, 20008,    -1, 20010, 20011, 20012, 20013, 20014, 20015,
+  -1  , 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, 20113, 20114, 20115,
+ 20100, 20101, 20102, 20103, 20104, 20105, 20106, 20107, 20108, 20109, 20110, 20111, 20112,    -1, 20114, 20115,
   //================== 405
- 20006,    -1,    -1, 20009,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1,    -1
-
+  //  0     1      2      3      4      5      6      7      8      9     10     11     12     13     14     15
+ 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;
@@ -141,8 +149,9 @@ void PrintMapping(){
 
         printf("\033[35m%3d(%2d)\033[0m|", detID, ch);
         
-      }else{
+      }else if( typeID == 3){ // MISC
 
+        printf("\033[33m%3d(%2d)\033[0m|", detID, ch);
 
       }
     }
@@ -215,7 +224,6 @@ 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") {

slope_intercept_cathode.txt

@@ -0,0 +1,23 @@
+Histogram Number	Slope	Intercept
+24	1	-2.89219e-10
+25	0.942098	-0.105169
+26	0.980862	-0.732032
+27	0.982975	-2.22704
+28	0.978815	-1.51477
+29	0.965245	-2.19515
+30	0.945384	-0.892599
+31	0.977408	-0.908592
+32	0.919546	3.25464
+33	0.972194	2.44956
+34	0.92852	5.44745
+35	0.947098	1.40531
+36	0.875491	-1.13145
+37	1.95496	-1735.58
+38	0.970862	2.86019
+40	0.91793	-3.80615
+41	0.913897	-2.12964
+42	0.954014	-0.760604
+43	0.993616	-1.40278
+45	0.926169	-21.2016
+46	1.00577	-2.14281
+47	0.943312	-1.26464

slope_intercept_results.txt

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

slope_intercept_results_anode.txt

@@ -0,0 +1,21 @@
+Histogram Number	Slope	Intercept
+1	1	-1.87356e-10
+2	0.964185	1.49989
+3	0.92638	-1.30621
+4	0.905569	1.00834
+5	0.901182	0.470903
+7	0.942785	1.08887
+8	0.878904	-0.0107433
+10	0.903343	8.38332
+11	0.914227	6.56108
+12	0.961008	23.0982
+13	0.920976	5.22104
+14	0.936584	31.5073
+15	0.959044	5.43267
+16	0.95263	-0.404053
+17	0.90953	4.82833
+18	0.940277	10.3629
+20	1.00683	4.76371
+21	0.968342	-43.9496
+22	0.892882	-32.0742
+23	0.933615	1.10704

sx3_BackGains.txt

@@ -0,0 +1,12 @@
+7 0 1.20298
+7 1 0.995493
+7 2 0.993613
+7 3 1.2514
+9 0 1.01574
+9 1 0.961032
+9 2 0.988379
+9 3 1.05832
+19 0 1.07936
+19 1 0.97626
+19 2 1.00078
+19 3 1.03335

sx3_GainMatchfront.txt

@@ -0,0 +1,8 @@
+3 0 1 1 0.852399 0.881228
+3 0 2 2 0.813845 0.975967
+3 0 3 3 0.859643 0.863715
+3 1 1 1 0.76728 0.942438
+3 1 2 2 0.780302 0.929008
+3 2 2 2 0.729082 1.02005
+3 3 2 2 0.759098 1.05376
+3 3 3 3 0.821183 0.952335