modified: Analyzer.C

modified:   ProcessRun.sh
	modified:   mapping.h

.vscode/settings.json

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

Analyzer.C

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

Analyzer.h

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

Analyzer1.C

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

Analyzer1.h

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

Armory/ClassDet.h

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

Armory/ClassPW.h

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

Armory/Mapper.cpp

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

PCGainMatch.C

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

PCGainMatch.h

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

PreAnalysis_C_ACLiC_dict.cxx_tmp_829464

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

mapping.h

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