modified: .gitignore

modified:   Analyzer.C

.gitignore

@@ -10,3 +10,4 @@ AnasenMS
 
 data/
 data_proton/
+root_data/

.vscode/settings.json

@@ -100,12 +100,7 @@
     "PCPulser_All_new.C": "cpp",
     "PosCal_2.C": "cpp",
     "AutoFit.C": "cpp",
-    "Fitting.C": "cpp",
+    "Fitting.C": "cpp"
-    "PCGainMatch.C": "cpp",
-    "Analyzer1.C": "cpp",
-    "FitHistogramsWithTSpectrum_Sequential_Improved.C": "cpp",
-    "PlotAndFitCentroids.C": "cpp",
-    "MatchAndPlotCentroids.C": "cpp"
   },
   "github-enterprise.uri": "https://fsunuc.physics.fsu.edu"
 }

Analysis.C

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

Analyzer.C

@@ -1,27 +1,23 @@
 #define Analyzer_cxx
 
 #include "Analyzer.h"
-#include "Armory/ClassSX3.h"
-#include "Armory/ClassPW.h"
-
 #include <TH2.h>
 #include <TStyle.h>
 #include <TCanvas.h>
 #include <TMath.h>
-#include "TVector3.h"
 
-#include <fstream>
-#include <iostream>
-#include <sstream>
-#include <map>
 #include <utility>
 #include <algorithm>
 
+#include "Armory/ClassSX3.h"
+#include "Armory/ClassPW.h"
+
+#include "TVector3.h"
+
 TH2F *hsx3IndexVE;
 TH2F *hqqqIndexVE;
 TH2F *hpcIndexVE;
 
-TH2F *hpcIndexVE_GM;
 TH2F *hsx3Coin;
 TH2F *hqqqCoin;
 TH2F *hpcCoin;
@@ -32,23 +28,15 @@ TH2F *hqqqVpcIndex;
 TH2F *hqqqVpcE;
 TH2F *hsx3VpcE;
 TH2F *hanVScatsum;
-TH2F *hanVScatsum_a[24];
-TH1F *hPC_E[48];
-TH1F *hAnodeMultiplicity;
-
 int padID = 0;
 
 SX3 sx3_contr;
 PW pw_contr;
 PW pwinstance;
 TVector3 hitPos;
-// TVector3 anodeIntersection;
-std::map<int, std::pair<double, double>> slopeInterceptMap;
-
 bool HitNonZero;
 
 TH1F *hZProj;
-TH1F *hPCZProj;
 
 void Analyzer::Begin(TTree * /*tree*/)
 {
@@ -58,10 +46,8 @@ void Analyzer::Begin(TTree * /*tree*/)
   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, 16000);
+  hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2 * 24, 0, 2 * 24, 400, 0, 4000);
   hpcIndexVE->SetNdivisions(-1204, "x");
-  hpcIndexVE_GM = new TH2F("hpcIndexVE_GM", "PC index vs Energy; PC index ; Energy", 2 * 24, 0, 2 * 24, 400, 0, 16000);
-  hpcIndexVE_GM->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);
@@ -84,50 +70,11 @@ void Analyzer::Begin(TTree * /*tree*/)
   hsx3VpcE->SetNdivisions(-612, "x");
   hsx3VpcE->SetNdivisions(-12, "y");
 
-  hZProj = new TH1F("hZProj", "Z Projection", 1200, -600, 600);
+  hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
-  hPCZProj = new TH1F("hPCZProj", "PC Z Projection", 1200, -600, 600);
 
-  hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400, 0, 16000, 400, 0, 20000);
+  hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400, 0, 10000, 400, 0, 16000);
-  hAnodeMultiplicity = new TH1F("hAnodeMultiplicity", "Number of Anodes/Event", 40, 0, 40);
-  // 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, 16000, 400, 0, 20000);
-  // }
-  for (int i = 0; i < 48; i++)
-  {
-    TString histName = Form("hCathode_%d", i);
-    TString histTitle = Form("Cathode_E_%d;", i);
-    hPC_E[i] = new TH1F(histName, histTitle, 3200, 0, 32000);
-  }
   sx3_contr.ConstructGeo();
   pw_contr.ConstructGeo();
-
-  std::ifstream inputFile("slope_intercept_results.txt");
-
-  if (inputFile.is_open())
-  {
-    std::string line;
-    int index;
-    double slope, intercept;
-    while (std::getline(inputFile, line))
-    {
-      std::stringstream ss(line);
-      ss >> index >> slope >> intercept;
-      // wires 37, 39, 44 have fit data that is incorrect or not present, they have thus been set to 1,0 (slope, intercept) for convenience
-      // wire 19 the 4th point was genereated using the slope of the line produced uising the other 3 points from the wire 1 vs wire 19 plot
-      if (index >= 0 && index <= 47)
-      {
-        slopeInterceptMap[index] = std::make_pair(slope, intercept);
-      }
-    }
-    inputFile.close();
-  }
-  else
-  {
-    std::cerr << "Error opening slope_intercept.txt" << std::endl;
-  }
 }
 
 Bool_t Analyzer::Process(Long64_t entry)
@@ -332,10 +279,29 @@ Bool_t Analyzer::Process(Long64_t entry)
   }
   // //======================= PC
 
-  // Calculate the crossover points and put them into an array
+  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 );
 
   pwinstance.ConstructGeo();
-  Coord Crossover[24][24][2];
+  Coord Crossover[24][24];
   TVector3 a, c, diff;
   double a2, ac, c2, adiff, cdiff, denom, alpha, beta;
   int index = 0;
@@ -344,11 +310,6 @@ Bool_t Analyzer::Process(Long64_t entry)
     a = pwinstance.An[i].first - pwinstance.An[i].second;
     for (int j = 0; j < pwinstance.Ca.size(); j++)
     {
-      // Ok so this method uses what is essentially the solution of 2 equations to find the point of intersection between the anode and cathode wires
-      // here a and c are the vectors of the anode and cathode wires respectively
-      // diff is the perpendicular vector between the anode and cathode wires
-      // The idea behind this is to then find the scalars alpha and beta that give a ratio between 0 and -1,
-
       c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
       diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
       a2 = a.Dot(a);
@@ -359,246 +320,82 @@ Bool_t Analyzer::Process(Long64_t entry)
       denom = a2 * c2 - ac * ac;
       alpha = (ac * cdiff - c2 * adiff) / denom;
       beta = (a2 * cdiff - ac * adiff) / denom;
-
+      Crossover[i][j].x = pwinstance.An[i].first.X() + alpha * a.X();
-      Crossover[i][j][0].x = pwinstance.An[i].first.X() + alpha * a.X();
+      Crossover[i][j].y = pwinstance.An[i].first.Y() + alpha * a.Y();
-      Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
+      Crossover[i][j].z = pwinstance.An[i].first.Z() + alpha * a.Z();
-      Crossover[i][j][0].z = pwinstance.An[i].first.Z() + alpha * a.Z();
+      if (i == 23)
-
-      // placeholder variable Crossover[i][j][2].x has nothing to do with the geometry of the crossover and is being used to store the alpha value-
-      //-so that it can be used to sort "good" hits later
-      Crossover[i][j][1].x = alpha;
-
-      if (i == 4)
       {
-        for (int k = -4; k < 3; k++)
+        if (abs(i - j) < 7 || abs(i - j) > 17)
         {
-          if ((i + 24 + k) % 24 == 23 - j) // the 23-j is used to accomodate for the fact that the order of the cathodes was reversed
+          if (alpha < 0 && alpha > -1)
           {
-            // if (alpha < 1 && alpha >= -1)
+            printf("Anode and cathode indices and coord : %d %d %f %f %f %f\n", i, j, pwinstance.Ca[j].first.X(), pwinstance.Ca[j].first.Y(), pwinstance.Ca[j].first.Z(), alpha);
-            {
+            printf("Crossover wires, points and alpha are : %f %f %f %f \n", Crossover[i][j].x, Crossover[i][j].y, Crossover[i][j].z, alpha);
-              Crossover[i][j][0].z = 9999999; // this is a placeholder value to indicate that the anode and cathode wires do not intersect
-              printf("Anode and cathode indices and coord : %d %d %f %f %f %f\n", i, j, pwinstance.Ca[j].first.X(), pwinstance.Ca[j].first.Y(), pwinstance.Ca[j].first.Z(), alpha);
-              // printf("Crossover wires, points and alpha are : %f %f %f %f \n", Crossover[i][j][1].x, Crossover[i][j][1].y, Crossover[i][j][1].z, Crossover[i][j][2].x /*this is alpha*/);
-            }
           }
         }
       }
     }
   }
 
-  // anodeIntersection.Clear();
+  if (E.size() >= 3)
-  for (int i = 0; i < pc.multi; i++)
   {
 
-    if (pc.e[i] > 100)
+    int aID = 0;
+    int cID = 0;
+
+    float aE = 0;
+    float cE = 0;
+    // 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++)
     {
-      hpcIndexVE->Fill(pc.index[i], pc.e[i]); // non gain matched energy
+      if (qqq.index[k] == 75 && pc.index[k] == 2 && pc.e[k] > 100)
-    }
-
-    // Gain Matching of PC wires
-    if (pc.index[i] >= 0 && pc.index[i] < 48)
-    {
-      // printf("index: %d, Old cathode energy: %d \n", pc.index[i],pc.e[i]);
-      auto it = slopeInterceptMap.find(pc.index[i]);
-      if (it != slopeInterceptMap.end())
       {
-        double slope = it->second.first;
-        double intercept = it->second.second;
-        // printf("slope: %f, intercept:%f\n" ,slope, intercept);
-        pc.e[i] = slope * pc.e[i] + intercept;
-        // printf("index: %d, New cathode energy: %d \n",pc.index[i], pc.e[i]);
-      }
-      hpcIndexVE_GM->Fill(pc.index[i], pc.e[i]);
-      // hPC_E[pc.index[i]]->Fill(pc.e[i]);
-    }
-  }
 
-  std::vector<std::pair<int, double>> anodeHits = {};
+        int multi_an = 0;
-  std::vector<std::pair<int, double>> cathodeHits = {};
+        for (int l = 0; l < E.size(); l++)
-  std::vector<std::pair<int, double>> corrcatMax = {};
-  std::vector<std::pair<int, double>> corrcatnextMax = {};
-  std::vector<std::pair<int, double>> commcat = {};
-  int aID = 0;
-  int cID = 0;
-  float aE = 0;
-  float cE = 0;
-  float aESum = 0;
-  float cESum = 0;
-  float aEMax = 0;
-  float cEMax = 0;
-  float aEnextMax = 0;
-  float cEnextMax = 0;
-  int aIDMax = 0;
-  int cIDMax = 0;
-  int aIDnextMax = 0;
-  int cIDnextMax = 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] > 100 /*&& pc.multi < 7*/)
-    {
-      // creating a vector of pairs of anode and cathode hits
-      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] - 24, 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]);
-      }
-    }
-  }
-  // sorting the anode and cathode hits in descending order of energy
-  std::sort(anodeHits.begin(), anodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
-            { return a.second > b.second; });
-  std::sort(cathodeHits.begin(), cathodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
-            { return a.second > b.second; });
-
-  if (anodeHits.size() >= 1 && cathodeHits.size() >= 1)
-  {
-
-    for (const auto &anode : anodeHits)
-    {
-      aID = anode.first;
-      aE = anode.second;
-      aESum += aE;
-      if (aE > aEMax)
-      {
-        aEMax = aE;
-        aIDMax = aID;
-      }
-      if (aE > aEnextMax && aE < aEMax)
-      {
-        aEnextMax = aE;
-        aIDnextMax = aID;
-      }
-    }
-
-    // std::cout << " Anode iD : " << aIDMax << " Energy : " << aEMax << std::endl;
-
-    // printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
-    for (const auto &cathode : cathodeHits)
-    {
-      cID = cathode.first;
-      cE = cathode.second;
-      // std::cout << "Cathode ID : " << cID << " Energy : " << cE << std::endl;
-      // if (cE > cEMax)
-      // {
-      //   cIDnextMax = cIDMax;
-      //   cEnextMax = cEMax;
-      //   cEMax = cE;
-      //   cIDMax = cID;
-      // }
-      // if (cE > cEnextMax && cE < cEMax)
-      // {
-      //   cEnextMax = cE;
-      //   cIDnextMax = cID;
-      // }
-
-      cESum += cE;
-      // This section of code is used to find the cathodes are correlated with the max and next max anodes, as well as to figure out if there are any common cathodes
-      // the anodes are correlated with the cathodes +/-2 from the anode number in the reverse order
-      for (int j = -3; j < 2; j++)
-      {
-        if ((aIDMax + 24 + j) % 24 == 23 - cID) /* the 23-cID is used to accomodate for the fact that the order of the cathodes was reversed relative top the physical geometry */
         {
-          corrcatMax.push_back(std::pair<int, double>(cID, cE));
+          if (E[l].first < 24)
-          // std::cout << " Cathode iD : " << cID << " Energy : " << cE << std::endl;
+          {
+            multi_an++;
+          }
+        }
+
+        if (multi_an >= 1)
+        {
+          for (int l = 0; l < E.size(); l++)
+          {
+            if (E[l].first < 24 && E[l].first != 19 && E[l].first != 12)
+            {
+              aE = E[l].second;
+            }
+            else if (E[l].first > 24)
+            {
+              cE = E[l].second;
+            }
+          }
         }
-        // if ((aIDnextMax + 24 + j) % 24 == cID)
-        // {
-        //   corrcatnextMax.push_back(std::pair<int, double>(cathode.first, cathode.second));
-        //   std::sort(corrcatMax.begin(), corrcatMax.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
-        //             { return a.second > b.second; });
-        // }
-        // for (int k = 0; k < 5; k++)
-        // {
-        //   if ((aIDMax + 24 + j) % 24 == (aIDnextMax + 24 + k) % 24)
-        //   {
-        //     commcat.push_back(std::pair<int, double>(cathode.first, cathode.second));
-        //   }
-        // }
       }
     }
-    // std::sort(corrcatMax.begin(), corrcatMax.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b ) { return a.second > b.second; });
+    hanVScatsum->Fill(aE, cE);
 
-    TVector3 anodeIntersection;
+    if (ID[0].first < 1)
-    // Implementing a method for PC reconstruction using a single Anode event
-    // if (anodeHits.size() == 1)
     {
-      for (const auto &corr : corrcatMax)
+      aID = pc.ch[ID[0].second];
-      {
+      cID = pc.ch[ID[1].second];
-        anodeIntersection += TVector3((corr.second) * 1.0 / cESum * Crossover[aIDMax][corr.first][0].x * 1.0, (corr.second) * 1.0 / cESum * Crossover[aIDMax][corr.first][0].y * 1.0,
+    }
-                                      (corr.second) * 1.0 / cESum * Crossover[aIDMax][corr.first][0].z * 1.0);
+    else
-        // std::cout << "Anode Intersection " << anodeIntersection.Z() << std::endl;
+    {
-      }
+      cID = pc.ch[ID[0].second];
+      aID = pc.ch[ID[1].second];
     }
-
-    // Filling the PC Z projection histogram
-    // std::cout << anodeIntersection.Z() << std::endl;
-    hPCZProj->Fill(anodeIntersection.Z());
-
-    // }
-
-    // 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(aEMax, 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(anodeHits.size()==1){
-    hAnodeMultiplicity->Fill(corrcatMax.size());
-    // }
 
     if (HitNonZero)
     {
@@ -693,23 +490,4 @@ void Analyzer::Terminate()
   canvas->cd(padID)->SetGrid(1);
   //  hZProj->Draw();
   hanVScatsum->Draw("colz");
-
-  TFile *outRoot = new TFile("Histograms.root", "RECREATE");
-
-  if (!outRoot->IsOpen())
-  {
-    std::cerr << "Error opening file for writing!" << std::endl;
-    return;
-  }
-
-  // Loop through histograms and write them to the ROOT file
-  for (int i = 0; i < 48; i++)
-  {
-    if (hPC_E[i] != nullptr)
-    {
-      hPC_E[i]->Write(); // Write histogram to file
-    }
-  }
-
-  outRoot->Close();
 }

Analyzer.h

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

Analyzer1.C

@@ -1,402 +0,0 @@
-#define Analyzer1_cxx
-
-#include "Analyzer1.h"
-#include <TH2.h>
-#include <TStyle.h>
-#include <TCanvas.h>
-#include <TMath.h>
-
-#include <utility>
-#include <algorithm>
-
-#include "Armory/ClassSX3.h"
-#include "Armory/ClassPW.h"
-
-#include "TVector3.h"
-
-TH2F * hsx3IndexVE;
-TH2F * hqqqIndexVE;
-TH2F * hpcIndexVE;
-
-TH2F * hsx3Coin;
-TH2F * hqqqCoin;
-TH2F * hpcCoin;
-
-TH2F * hqqqPolar;
-TH2F * hsx3VpcIndex;
-TH2F * hqqqVpcIndex;
-TH2F * hqqqVpcE;
-TH2F * hsx3VpcE;
-TH2F * hanVScatsum;
-int padID = 0;
-
-SX3 sx3_contr;
-PW pw_contr;
-TVector3 hitPos;
-bool HitNonZero;
-
-TH1F * hZProj;
-
-void Analyzer1::Begin(TTree * /*tree*/){
-  TString option = GetOption();
-
-  hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy",  24*12, 0,  24*12, 400, 0, 5000); hsx3IndexVE->SetNdivisions( -612, "x");
-  hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
-  hpcIndexVE = new TH2F("hpcIndexVE",   "PC index vs Energy; PC index ; Energy",     2*24, 0,   2*24, 400, 0, 4000); hpcIndexVE->SetNdivisions( -1204, "x");
-
-
-  hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident",  24*12, 0,  24*12,  24*12, 0, 24*12);
-  hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4*2*16, 0, 4*2*16, 4*2*16, 0, 4*2*16);
-  hpcCoin  = new TH2F("hpcCoin",  "PC Coincident",    2*24, 0,   2*24,   2*24, 0, 2*24);
-
-  hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16*4, -TMath::Pi(), TMath::Pi(),16, 10, 50);
-
-  hsx3VpcIndex = new TH2F("hsx3Vpcindex", "sx3 vs pc; sx3 index; pc index",  24*12, 0, 24*12, 48, 0, 48);
-  hsx3VpcIndex->SetNdivisions( -612, "x");
-  hsx3VpcIndex->SetNdivisions( -12, "y");
-  hqqqVpcIndex = new TH2F("hqqqVpcindex", "qqq vs pc; qqq index; pc index",  4*2*16, 0, 4*2*16, 48, 0, 48);
-  hqqqVpcIndex->SetNdivisions( -612, "x");
-  hqqqVpcIndex->SetNdivisions( -12, "y");
-
-  hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 400, 0, 5000);
-  hqqqVpcE->SetNdivisions( -612, "x");
-  hqqqVpcE->SetNdivisions( -12, "y");
-  
-  hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000);
-  hsx3VpcE->SetNdivisions( -612, "x");
-  hsx3VpcE->SetNdivisions( -12, "y");
-
-  hZProj = new TH1F("hZProj", "Z Projection", 1200, -600, 600);
-
-  hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000);
-
-  sx3_contr.ConstructGeo();
-  pw_contr.ConstructGeo();
-
-}
-
-
-
-
-Bool_t Analyzer1::Process(Long64_t entry){
-
-  // if ( entry > 100 ) return kTRUE;
-
-  hitPos.Clear();
-  HitNonZero = false;
-
-  // if( entry > 1) return kTRUE;
-  // printf("################### ev : %llu \n", entry);
-
-  b_sx3Multi->GetEntry(entry);
-  b_sx3ID->GetEntry(entry);
-  b_sx3Ch->GetEntry(entry);
-  b_sx3E->GetEntry(entry);
-  b_sx3T->GetEntry(entry);
-  b_qqqMulti->GetEntry(entry);
-  b_qqqID->GetEntry(entry);
-  b_qqqCh->GetEntry(entry);
-  b_qqqE->GetEntry(entry);
-  b_qqqT->GetEntry(entry);
-  b_pcMulti->GetEntry(entry);
-  b_pcID->GetEntry(entry);
-  b_pcCh->GetEntry(entry);
-  b_pcE->GetEntry(entry);
-  b_pcT->GetEntry(entry);
-
-  sx3.CalIndex();
-  qqq.CalIndex();
-  pc.CalIndex();
-
-  // sx3.Print();
-
-  //########################################################### Raw data
-  // //======================= SX3
-
-  std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
-  for( int i = 0; i < sx3.multi; i ++){
-    ID.push_back(std::pair<int, int>(sx3.id[i], i));
-
-    hsx3IndexVE->Fill( sx3.index[i], sx3.e[i] );
-
-    for( int j = i+1; j < sx3.multi; j++){
-      hsx3Coin->Fill( sx3.index[i], sx3.index[j]);
-    }
-
-    for( int j = 0; j < pc.multi; j++){
-      hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
-      //  if( sx3.ch[index] > 8 ){
-      //    hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
-      //   }
-    }
-  }
-
-
-  if( ID.size() > 0 ){
-    std::sort(ID.begin(), ID.end(),  [](const std::pair<int, int> & a, const std::pair<int, int> & b) {
-      return a.first < b.first;
-    } );
-    // printf("##############################\n");
-    // for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second );
-
-    std::vector<std::pair<int, int>> sx3ID; 
-    sx3ID.push_back(ID[0]);
-    bool found = false;
-    for( size_t i = 1; i < ID.size(); i++){
-      if( ID[i].first == sx3ID.back().first) {
-        sx3ID.push_back(ID[i]);
-        if( sx3ID.size() >= 3) {
-          found = true;
-        } 
-      }else{
-        if( !found ){
-          sx3ID.clear();
-          sx3ID.push_back(ID[i]);
-        }
-      }
-    }
-
-    // printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
-
-    if( found ){
-      int sx3ChUp, sx3ChDn, sx3ChBk;
-      float sx3EUp, sx3EDn;
-      // printf("------ sx3 ID : %d, multi: %zu\n", sx3ID[0].first, sx3ID.size());
-      for( size_t i = 0; i < sx3ID.size(); i++ ){
-        int index = sx3ID[i].second;
-        // printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]);
-
-
-        if( sx3.ch[index] < 8 ){
-          if( sx3.ch[index] % 2 == 0) {
-            sx3ChDn = sx3.ch[index];
-            sx3EDn = sx3.e[index];
-          }else{
-            sx3ChUp = sx3.ch[index];
-            sx3EUp = sx3.e[index];
-          }
-        }else{
-          sx3ChBk = sx3.ch[index];
-        }
-         for( int j = 0; j < pc.multi; j++){
-      // hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
-            if( sx3.ch[index] > 8 ){
-              hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
-          //  hpcIndexVE->Fill( pc.index[i], pc.e[i] );
-            }
-          }
-      }
-
-      sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
-      hitPos = sx3_contr.GetHitPos();
-      HitNonZero = true;
-      // hitPos.Print();
-    }
-
-  }
-  
-
-  // //======================= QQQ
-  for( int i = 0; i < qqq.multi; i ++){
-    // for( int j = 0; j < pc.multi; j++){
-        // if(pc.index[j]==4){
-          hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] );
-        // }
-    // }
-    for( int j = 0; j < qqq.multi; j++){
-      if ( j == i ) continue;
-      hqqqCoin->Fill( qqq.index[i], qqq.index[j]);
-    }
-
-
-    for( int j = i + 1; j < qqq.multi; j++){
-      for( int k = 0; k < pc.multi; k++){
-        if(pc.index[k]<24 && pc.e[k]>50 ){
-          hqqqVpcE->Fill( qqq.e[i], pc.e[k] );
-          //  hpcIndexVE->Fill( pc.index[i], pc.e[i] );
-                hqqqVpcIndex->Fill( qqq.index[i], pc.index[j] );
-
-        }
-      // }
-      }
-      // if( qqq.used[i] == true ) continue;
-      
-      //if( qqq.id[i] == qqq.id[j] && (16 - qqq.ch[i]) * (16 - qqq.ch[j]) < 0  ){ // must be same detector and wedge and ring
-      if( qqq.id[i] == qqq.id[j]  ){ // must be same detector 
-        
-        int chWedge = -1;
-        int chRing  =  -1;
-        if( qqq.ch[i] < qqq.ch[j]){
-          chRing = qqq.ch[j] - 16;
-          chWedge = qqq.ch[i];
-        }else{
-          chRing = qqq.ch[i];
-          chWedge = qqq.ch[j] - 16;
-        }
-
-        // printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing);
-
-        double theta = -TMath::Pi()/2 + 2*TMath::Pi()/16/4.*(qqq.id[i]*16 + chWedge +0.5);
-        double rho   = 10.+40./16.*(chRing+0.5);
-        // if(qqq.e[i]>50){
-        hqqqPolar->Fill( theta, rho);
-        // }
-        // qqq.used[i] = true;
-        // qqq.used[j] = true;
-
-        if( !HitNonZero ){
-        double x = rho * TMath::Cos(theta);
-        double y = rho * TMath::Sin(theta);
-        hitPos.SetXYZ(x, y, 23 + 75 + 30);
-        HitNonZero = true;
-        }
-      }
-    }
-  
-   
-  }
-  // //======================= PC
-
-  ID.clear();
-  int counter=0;
-  std::vector<std::pair<int, double>> E; 
-  E.clear();
-  for( int i = 0; i < pc.multi; i ++){
-
-    if( pc.e[i] > 100 ) ID.push_back(std::pair<int, int>(pc.id[i], i));
-    if( pc.e[i] > 100 ) E.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
-
-    hpcIndexVE->Fill( pc.index[i], pc.e[i] );
-
-    for( int j = i+1; j < pc.multi; j++){
-      hpcCoin->Fill( pc.index[i], pc.index[j]);
-    
-    }
-
-  }
-  //  for( size_t i = 0; i < E.size(); i++) printf("%zu | %d %d \n", i, E[i].first, E[i].second ); 
-
-  if( E.size()>=3 ){
-
-    int aID = 0;
-    int cID = 0;
-
-    float aE = 0;
-    float cE = 0;
-    bool multi_an =false;
-    // if( ID[0].first < 1 ) {
-    //   aID = pc.ch[ID[0].second];
-    //   cID = pc.ch[ID[1].second];
-    // }else{
-    //   cID = pc.ch[ID[0].second];
-    //   aID = pc.ch[ID[1].second];
-    // }
-    // printf("anode= %d, cathode = %d\n", aID, cID);
-
-  // for( int k = 0; k < qqq.multi; k++){
-  //   if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){
-      for(int l=0;l<E.size();l++){
-        if(E[l].first<24 ){
-          if(!multi_an){
-            aE = E[l].second;
-          }
-          multi_an=true;
-        }
-        else if (E[l].first>=24){
-          cE = E[l].second + cE;
-        }
-      }
-    // }
-  // }
-    hanVScatsum->Fill(aE,cE);
-      
-    if( ID[0].first < 1 ) {
-      aID = pc.ch[ID[0].second];
-      cID = pc.ch[ID[1].second];
-    }else{
-      cID = pc.ch[ID[0].second];
-      aID = pc.ch[ID[1].second];
-    }
-
-    if( HitNonZero){
-      pw_contr.CalTrack( hitPos, aID, cID);
-      hZProj->Fill(pw_contr.GetZ0());
-    }
-
-  
-  }
-  
-
-
-  //########################################################### Track constrcution
-
-
-  //############################## DO THE KINEMATICS
-
-
-  return kTRUE;
-}
-
-void Analyzer1::Terminate(){
-
-  gStyle->SetOptStat("neiou");
-  TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); 
-  canvas->Divide(3,3);
-
-  //hsx3VpcIndex->Draw("colz");
-
-  //=============================================== pad-1
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hsx3IndexVE->Draw("colz");
-
-  //=============================================== pad-2
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hqqqIndexVE->Draw("colz");
-
-  //=============================================== pad-3
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hpcIndexVE->Draw("colz");
-
-  //=============================================== pad-4
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hsx3Coin->Draw("colz");
-
-  //=============================================== pad-5
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  canvas->cd(padID)->SetLogz(true);
-
-  hqqqCoin->Draw("colz");
-
-  //=============================================== pad-6
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hpcCoin->Draw("colz");
-
-  //=============================================== pad-7
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  // hsx3VpcIndex ->Draw("colz"); 
-  hsx3VpcE->Draw("colz") ;
-
-  //=============================================== pad-8
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  // hqqqVpcIndex ->Draw("colz");  
-
-  hqqqVpcE ->Draw("colz");
-  //=============================================== pad-9
-  padID ++; 
-
-  // canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
-  // hqqqPolar->Draw("same colz pol");
-
- canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-//  hZProj->Draw();
-  hanVScatsum->Draw("colz");
-
-}

Armory/ClassPW.h

@@ -163,10 +163,6 @@ inline void PW::ConstructGeo()
                      -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
-  std::rotate(Ca.begin(), Ca.begin() + 3, Ca.end());
 
   dAngle = wireShift * TMath::TwoPi() / nWire;
   anodeLength = TMath::Sqrt(zLen * zLen + TMath::Power(2 * radiusA * TMath::Sin(dAngle / 2), 2));

FitHistogramsWithTSpectrum_Sequential_Improved.C

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

MatchAndPlotCentroids.C

@@ -1,133 +0,0 @@
-#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

@@ -1,454 +0,0 @@
-#define PCGainMatch_cxx
-
-#include "PCGainMatch.h"
-#include <TH2.h>
-#include <TStyle.h>
-#include <TCanvas.h>
-#include <TMath.h>
-#include <TCutG.h>
-#include <utility>
-#include <algorithm>
-#include "Armory/ClassSX3.h"
-#include "Armory/ClassPW.h"
-
-#include "TVector3.h"
-
-TH2F * hsx3IndexVE;
-TH2F * hqqqIndexVE;
-TH2F * hpcIndexVE;
-
-TH2F * hsx3Coin;
-TH2F * hqqqCoin;
-TH2F * hpcCoin;
-
-TH2F * hqqqPolar;
-TH2F * hsx3VpcIndex;
-TH2F * hqqqVpcIndex;
-TH2F * hqqqVpcE;
-TH2F * hsx3VpcE;
-TH2F * hanVScatsum;
-TH2F * hanVScatsum_a[24];
-TH2F * hanVScatsum_hcut;
-TH2F * hanVScatsum_lcut;
-TH2F * hAnodeHits;
-TH1F * hAnodeMultiplicity;
-
-
-int padID = 0;
-
-SX3 sx3_contr;
-PW pw_contr;
-TVector3 hitPos;
-bool HitNonZero;
-
-TH1F * hZProj;
-TCutG *AnCatSum_high;
-TCutG *AnCatSum_low;
-TCutG *PCCoinc_cut1;
-TCutG *PCCoinc_cut2;
-bool inCuth;
-bool inCutl;
-bool inPCCut;
-
-void PCGainMatch::Begin(TTree * /*tree*/){
-  TString option = GetOption();
-
-  hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy",  24*12, 0,  24*12, 400, 0, 5000); hsx3IndexVE->SetNdivisions( -612, "x");
-  hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
-  hpcIndexVE = new TH2F("hpcIndexVE",   "PC index vs Energy; PC index ; Energy",     2*24, 0,   2*24, 800, 0, 16000); hpcIndexVE->SetNdivisions( -1204, "x");
-
-
-  hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident",  24*12, 0,  24*12,  24*12, 0, 24*12);
-  hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4*2*16, 0, 4*2*16, 4*2*16, 0, 4*2*16);
-  hpcCoin  = new TH2F("hpcCoin",  "PC Coincident",    2*24, 0,   2*24,   2*24, 0, 2*24);
-
-  hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16*4, -TMath::Pi(), TMath::Pi(),16, 10, 50);
-
-  hsx3VpcIndex = new TH2F("hsx3Vpcindex", "sx3 vs pc; sx3 index; pc index",  24*12, 0, 24*12, 48, 0, 48);
-  hsx3VpcIndex->SetNdivisions( -612, "x");
-  hsx3VpcIndex->SetNdivisions( -12, "y");
-  hqqqVpcIndex = new TH2F("hqqqVpcindex", "qqq vs pc; qqq index; pc index",  4*2*16, 0, 4*2*16, 48, 0, 48);
-  hqqqVpcIndex->SetNdivisions( -612, "x");
-  hqqqVpcIndex->SetNdivisions( -12, "y");
-
-  hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 400, 0, 5000);
-  hqqqVpcE->SetNdivisions( -612, "x");
-  hqqqVpcE->SetNdivisions( -12, "y");
-  
-  hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000);
-  hsx3VpcE->SetNdivisions( -612, "x");
-  hsx3VpcE->SetNdivisions( -12, "y");
-
-  hZProj = new TH1F("hZProj", "Nos of anodes", 20, 0, 19);
-  hAnodeHits = new TH2F("hAnodeHits", "Anode vs Anode Energy, Anode ID; Anode E", 24,0 , 23, 400, 0 , 20000);
-  hAnodeMultiplicity = new TH1F("hAnodeMultiplicity", "Number of Anodes/Event", 40, 0, 40); 
-  hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 800, 0 , 16000);
-  for (int i = 0; i < 24; i++) {
-        TString histName = Form("hAnodeVsCathode_%d", i);
-        TString histTitle = Form("Anode %d vs Cathode Sum; Anode E; Cathode Sum E", i);
-        hanVScatsum_a[i] = new TH2F(histName, histTitle, 400, 0, 10000, 400, 0, 16000);
-    }
-  hanVScatsum_lcut = new TH2F("hanVScatsum_LCUT", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 16000, 400, 0 , 16000);
-  hanVScatsum_hcut = new TH2F("hanVScatsum_HCUT", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 16000, 400, 0 , 16000);
-
-  sx3_contr.ConstructGeo();
-  pw_contr.ConstructGeo();
-  // TFile *f1 = new TFile("AnCatSum_high.root");
-  // TFile *f2 = new TFile("AnCatSum_low.root");
-  // TFile *f3 = new TFile("PCCoinc_cut1.root");
-  // TFile *f4 = new TFile("PCCoinc_cut2.root");
-  // AnCatSum_high= (TCutG*)f1->Get("AnCatSum_high");
-  // AnCatSum_low= (TCutG*)f2->Get("AnCatSum_low");
-  // PCCoinc_cut1= (TCutG*)f3->Get("PCCoinc_cut1");
-  // PCCoinc_cut2= (TCutG*)f4->Get("PCCoinc_cut2");
-}
-
-
-
-
-Bool_t PCGainMatch::Process(Long64_t entry){
-// if (entry % 1000000 == 0) {
-//     std::cout << "Processing entry: " << entry << std::endl;
-// }
-  // if ( entry > 100 ) return kTRUE;
-
-  hitPos.Clear();
-  HitNonZero = false;
-
-  // if( entry > 1) return kTRUE;
-  // printf("################### ev : %llu \n", entry);
-
-  b_sx3Multi->GetEntry(entry);
-  b_sx3ID->GetEntry(entry);
-  b_sx3Ch->GetEntry(entry);
-  b_sx3E->GetEntry(entry);
-  b_sx3T->GetEntry(entry);
-  b_qqqMulti->GetEntry(entry);
-  b_qqqID->GetEntry(entry);
-  b_qqqCh->GetEntry(entry);
-  b_qqqE->GetEntry(entry);
-  b_qqqT->GetEntry(entry);
-  b_pcMulti->GetEntry(entry);
-  b_pcID->GetEntry(entry);
-  b_pcCh->GetEntry(entry);
-  b_pcE->GetEntry(entry);
-  b_pcT->GetEntry(entry);
-
-  sx3.CalIndex();
-  qqq.CalIndex();
-  pc.CalIndex();
-
-  // sx3.Print();
-
-  //########################################################### Raw data
-
-  // //======================= PC
-
-  std::vector<std::pair<int, double>> anodeHits={};
-  std::vector<std::pair<int, double>> cathodeHits={};  
-  int aID = 0;
-  int cID = 0;
-  float aE = 0;
-  float cE = 0;
-
-// Define the excluded SX3 and QQQ channels
-  // std::unordered_set<int> excludeSX3 = {34, 35, 36, 37, 61, 62, 67, 73, 74, 75, 76, 77, 78, 79, 80, 93, 97, 100, 103, 108, 109, 110, 111, 112};
-  // std::unordered_set<int> excludeQQQ = {0, 17, 109, 110, 111, 112, 113, 119, 127, 128};
-  // inCuth=false;
-  // inCutl=false;
-  // inPCCut=false;
-  for( int i = 0; i < pc.multi; i ++){
-    
-    if(pc.e[i]>50 && pc.multi<7){ 
-
-      float aESum = 0;
-      float cESum = 0;
-      if (pc.index[i] < 24 ) {
-            anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
-        } else if (pc.index[i] >= 24) {
-            cathodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
-        }
-
-      for(int j=i+1;j<pc.multi;j++){
-        // if(PCCoinc_cut1->IsInside(pc.index[i], pc.index[j]) || PCCoinc_cut2->IsInside(pc.index[i], pc.index[j])){
-        //   // hpcCoin->Fill(pc.index[i], pc.index[j]);
-        //   inPCCut = true;
-        // }
-        hpcCoin->Fill(pc.index[i], pc.index[j]);
-      } 
-      if (anodeHits.size()==1 && cathodeHits.size() >= 1) {
-          
-          for (const auto& anode : anodeHits) {
-          // for(int l=0; l<sx3.multi; l++){
-          //   if (sx3.index[l]==80){
-
-              aID = anode.first;
-              aE = anode.second;
-              aESum += aE;
-              printf("aID : %d, aE : %f\n", aID, aE);
-          }
-
-            printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
-            for (const auto& cathode : cathodeHits) {
-              cID = cathode.first;
-              cE = cathode.second;
-              // if(cE>cEMax){
-              //   cEMax = cE;
-              //   cIDMax = cID;
-              // }
-              // if(cE>cEnextMax && cE<cEMax){
-              //   cEnextMax = cE;
-              //   cIDnextMax = cID;
-              // }
-
-              cESum += cE;
-              }
-          // }
-          
-          // inCuth = false;
-          // inCutl = false;
-          // inPCCut = false;
-          // for(int j=i+1;j<pc.multi;j++){
-          //   if(PCCoinc_cut1->IsInside(pc.index[i], pc.index[j]) || PCCoinc_cut2->IsInside(pc.index[i], pc.index[j])){
-          //     // hpcCoin->Fill(pc.index[i], pc.index[j]);
-          //     inPCCut = true;
-          //   }
-          //   hpcCoin->Fill(pc.index[i], pc.index[j]);
-          // } 
-
-          // Check if the accumulated energies are within the defined ranges
-          // if (AnCatSum_high && AnCatSum_high->IsInside(aESum, cESum)) {
-          //     inCuth = true;
-          // }
-          // if (AnCatSum_low && AnCatSum_low->IsInside(aESum, cESum)) {
-          //     inCutl = true;
-          // }
-          
-          // Fill histograms based on the cut conditions
-          // if (inCuth && inPCCut) {
-          //     hanVScatsum_hcut->Fill(aESum, cESum);
-          // }
-          // if (inCutl && inPCCut) {
-          //     hanVScatsum_lcut->Fill(aESum, cESum);
-          // }               
-            // for(auto anode : anodeHits){
-         
-              // float aE = anode.second;
-              // aESum += aE;
-              // if(inPCCut){
-            hanVScatsum->Fill(aESum, cESum);
-            // }
-            if (aID < 24 && aE > 50) {
-             hanVScatsum_a[aID]->Fill(aE, cESum);
-            }
-
-          
-            // }
-          // Fill histograms for the `pc` data
-          hpcIndexVE->Fill(pc.index[i], pc.e[i]);
-          // if(inPCCut){
-            hAnodeMultiplicity->Fill(anodeHits.size());
-          // }
-      }
-
-    }
-  }    
-  
-
-
-  // //======================= SX3
-
-  std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
-  for( int i = 0; i < sx3.multi; i ++){
-    if(sx3.e[i]>50){
-      ID.push_back(std::pair<int, int>(sx3.id[i], i));
-      hsx3IndexVE->Fill( sx3.index[i], sx3.e[i] );
-
-      for( int j = i+1; j < sx3.multi; j++){
-        hsx3Coin->Fill( sx3.index[i], sx3.index[j]);
-      }
-
-      for( int j = 0; j < pc.multi; j++){
-        hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
-        //  if( sx3.ch[index] > 8 ){
-        //    hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
-        //   }
-      }
-    } 
-  }
-
-
-  if( ID.size() > 0 ){
-    std::sort(ID.begin(), ID.end(),  [](const std::pair<int, int> & a, const std::pair<int, int> & b) {
-      return a.first < b.first;
-    } );
-    // printf("##############################\n");
-    // for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second );
-
-    std::vector<std::pair<int, int>> sx3ID; 
-    sx3ID.push_back(ID[0]);
-    bool found = false;
-    for( size_t i = 1; i < ID.size(); i++){
-      if( ID[i].first == sx3ID.back().first) {
-        sx3ID.push_back(ID[i]);
-        if( sx3ID.size() >= 3) {
-          found = true;
-        } 
-      }else{
-        if( !found ){
-          sx3ID.clear();
-          sx3ID.push_back(ID[i]);
-        }
-      }
-    }
-
-    // printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
-
-    if( found ){
-      int sx3ChUp, sx3ChDn, sx3ChBk;
-      float sx3EUp, sx3EDn;
-      // printf("------ sx3 ID : %d, multi: %zu\n", sx3ID[0].first, sx3ID.size());
-      for( size_t i = 0; i < sx3ID.size(); i++ ){
-        int index = sx3ID[i].second;
-        // printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]);
-
-
-        if( sx3.ch[index] < 8 ){
-          if( sx3.ch[index] % 2 == 0) {
-            sx3ChDn = sx3.ch[index];
-            sx3EDn = sx3.e[index];
-          }else{
-            sx3ChUp = sx3.ch[index];
-            sx3EUp = sx3.e[index];
-          }
-        }else{
-          sx3ChBk = sx3.ch[index];
-        }
-         for( int j = 0; j < pc.multi; j++){
-
-      // hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
-            if( sx3.ch[index] > 8 && pc.index[j]<24 && pc.e[j]>50 ){
-              hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
-          //  hpcIndexVE->Fill( pc.index[i], pc.e[i] );
-            }
-          }
-      }
-
-      sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
-      hitPos = sx3_contr.GetHitPos();
-      HitNonZero = true;
-      // hitPos.Print();
-    }
-
-  }
-  
-
-  // //======================= QQQ
-  for( int i = 0; i < qqq.multi; i ++){
-
-      // for( int j = 0; j < pc.multi; j++){
-        if(qqq.e[i]>50 ){
-          hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] );
-          for( int j = 0; j < qqq.multi; j++){
-            if ( j == i ) continue;
-            hqqqCoin->Fill( qqq.index[i], qqq.index[j]);
-          }
-
-
-          for( int j = i + 1; j < qqq.multi; j++){
-            for( int k = 0; k < pc.multi; k++){
-              // if(qqq.e[i>50]){
-                hqqqVpcE->Fill( qqq.e[i], pc.e[k] );
-                hqqqVpcIndex->Fill( qqq.index[i], pc.index[j] );
-              }
-            // }
-          }
-        }  
-      // }
-  }
-
-    // hanVScatsum->Fill(aE,cE);
-    
-
-    if( HitNonZero){
-      pw_contr.CalTrack( hitPos, aID, cID);
-      hZProj->Fill(pw_contr.GetZ0());
-    }
-
-  
-  
-  
-
-
-  //########################################################### Track constrcution
-
-
-  //############################## DO THE KINEMATICS
-
-
-  return kTRUE;
-}
-
-void PCGainMatch::Terminate(){
-
-  gStyle->SetOptStat("neiou");
-  TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); 
-  canvas->Divide(3,3);
-
-  //hsx3VpcIndex->Draw("colz");
-
-  //=============================================== pad-1
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hsx3IndexVE->Draw("colz");
-
-  //=============================================== pad-2
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hqqqIndexVE->Draw("colz");
-
-  //=============================================== pad-3
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hpcIndexVE->Draw("colz");
-
-  //=============================================== pad-4
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hsx3Coin->Draw("colz");
-
-  //=============================================== pad-5
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  canvas->cd(padID)->SetLogz(true);
-
-  hqqqCoin->Draw("colz");
-
-  //=============================================== pad-6
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  hpcCoin->Draw("colz");
-
-  //=============================================== pad-7
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  // hsx3VpcIndex ->Draw("colz"); 
-  hsx3VpcE->Draw("colz") ;
-
-  //=============================================== pad-8
-  padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-
-  // hqqqVpcIndex ->Draw("colz");  
-
-  hqqqVpcE ->Draw("colz");
-  //=============================================== pad-9
-  padID ++; 
-
-  // canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
-  // hqqqPolar->Draw("same colz pol");
-
- canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
-  hanVScatsum->Draw("colz");
-  // hAnodeHits->Draw("colz");
-  // hAnodeMultiplicity->Draw();
-}

PreAnalysis_C_ACLiC_dict.cxx_tmp_829464

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

ProcessRun.sh

@@ -1,8 +1,7 @@
 #!/bin/bash
 
-if [ "$#" -ne 3 ]; then
+if [ "$#" -ne 2 ]; then
-    echo "Usage: $0 runID timeWindow_ns option"
+    echo "Usage: $0 runID timeWindow_ns"
-    echo "option: 0 - process raw data, 1 - process mapped data"
     echo "Exiting..."
     exit 1
 fi
@@ -10,9 +9,7 @@ fi
 runID=$1
 timeWindow=$2
 
-option=$3
+rawFolder=/home/tandem/Desktop/analysis/data
-
-rawFolder=/home/tandem/data1/2024_09_17Fap/data
 rootFolder=../root_data
 
 if [ $option -eq 0 ]; then
@@ -21,7 +18,7 @@ if [ $option -eq 0 ]; then
 
     fileList=`\ls -1 ${rawFolder}/*Run_${runID}_*.fsu`
 
-    ./EventBuilder ${timeWindow} 0 0 100000000 ${fileList}
+    ./EventBuilder ${timeWindow} 0 0 100000000 ${fileList} 
 
     outFile=${rawFolder}/*${runID}*${timeWindow}.root
 
@@ -29,5 +26,4 @@ if [ $option -eq 0 ]; then
 
     ./Mapper ${rootFolder}/*${runID}*${timeWindow}.root
 fi
-
+root "processRun.C(\"${rootFolder}/Run_${runID}_mapped.root\")"
-root "processRun.C(\"${rootFolder}/ProtonRun_${runID}_mapped.root\")"

TrackRecon.C

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

TrackRecon.h

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

centroids.txt

@@ -1,97 +0,0 @@
-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

@@ -1,89 +0,0 @@
-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

gainmatch.C

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

gainmatch.h

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

mapping.h

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

slope_intercept_cathode.txt

@@ -1,23 +0,0 @@
-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

@@ -1,49 +0,0 @@
-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

@@ -1,21 +0,0 @@
-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