/******************************************
 * compile with this:
 g++ -o JakeAnalyzer JakeAnalyzer.cpp `root-config --cflags --libs`
 */
#include <iostream>
#include <cmath>  // for std::isfinite
#include <TMath.h>  // for TMath::QuietNaN()
#include <fstream>  // For file handling
#include <sstream>  // For string stream parsing
//some root headers
#include <TTree.h>
#include <TH2F.h>
#include <TFile.h>
#include <TROOT.h>
#include <TString.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TChain.h>
#include <TSystem.h>
#include <TString.h>
#include <TCutG.h>
#include <TKey.h>
#include <regex>
//progress bar variables
const int barWidth = 50;
const char spinner[] = {'|', '/', '-', '\\'};
int spinIdx = 0;

// Detector setup
const Int_t nClovers = 11;
const Int_t nCrystals = (nClovers * 4);
const Int_t NGAGG = 64;
const Int_t nZeroDeg  = 2;
const Int_t nLaBr3 = 16;
const Int_t nBEGe = 2; // we recorded both of the channels off the bege
const Int_t nCrystalPairs = 66;
const Int_t nLabrPairs = 120;


//calibration variables// Global calibration arrays
Double_t ge_slope[nCrystals + nBEGe];
Double_t ge_intercept[nCrystals + nBEGe];
Double_t geTimeGates[nCrystalPairs];
Double_t labr_slope[nLaBr3];
Double_t labr_intercept[nLaBr3];


// PID variables
double traceSum[NGAGG];
double peak[NGAGG];
double tail[NGAGG];
int PIDxSize = 4000;
int PIDySize = 150000;
int PIDnBins = 1200;
Int_t coinCounter = 0;
int gaggID = 0;
double sumAB = 0;
double pidX = 0;
double pidY = 0;

//Comp Sup checker
bool bgoPresent = false;
Long64_t timeDiff_glob = 0;

//Ge Timing pair map file
const char* geTimingPairFile= "hatsuneMiku.dat";
//Labr Timing Pair Map File
const char* labrTimingPairFile= "hatsuneMiku2.dat";

//Global to hold the pairNum
Int_t pairNum = 0;

//global to hold the alphaCut check
bool alphaCheck = false;



Int_t makeGeTimingPairs(const char* filename, Int_t det1, Int_t det2) {
    std::ifstream infile(filename);
    if (!infile.is_open()) {
        std::cerr << "Error: Unable to open pairing file " << filename << std::endl;
        return -1;
    }

    int file_det1, file_det2, pairID;
    while (infile >> file_det1 >> file_det2 >> pairID) {
        if ((file_det1 == det1 && file_det2 == det2) ||
            (file_det1 == det2 && file_det2 == det1)) {
                //std::cout<<"det1: "<<file_det1<<" det2: "<<file_det2<<" pairID: "<<pairID<<std::endl;
                return pairID;
                
            }
        }

    std::cerr << "Pair (" << det1 << ", " << det2 << ") not found in " << filename << std::endl;
    return -1;
}

Int_t makeLabrTimingPairs(const char* filename, Int_t det1, Int_t det2) {
    std::ifstream infile(filename);
    if (!infile.is_open()) {
        std::cerr << "Error: Unable to open pairing file " << filename << std::endl;
        return -1;
    }

    int file_det1, file_det2, pairID;
    while (infile >> file_det1 >> file_det2 >> pairID) {
        if ((file_det1 == det1 && file_det2 == det2) ||
            (file_det1 == det2 && file_det2 == det1)) {
                //std::cout<<"det1: "<<file_det1<<" det2: "<<file_det2<<" pairID: "<<pairID<<std::endl;
                return pairID;
                
            }
        }

    std::cerr << "Pair (" << det1 << ", " << det2 << ") not found in " << filename << std::endl;
    return -1;
}

//Funtion to fill the PIDs
void makePID(TTree* tree, int kMaxMulti, int NGAGG, Int_t multi, Int_t* detNum, UInt_t (*qdc)[8], TH2F** hPID, bool trinCoin) {
    coinCounter = 0;
    for (int i = 0; i < NGAGG; ++i){
        for (int j = 0; j < 2; ++j){
            traceSum[i] = peak[i] = tail[i] = TMath::QuietNaN();
        }
    }
    //now enforcing coincidence between A and B channels
    for (int i = 0; i < multi; i++){
        if (i >= kMaxMulti || detNum[i] < 200 || detNum[i] > 399) continue;
        for (int j = 0; j < multi; j++){
            if (detNum[i]<=299){
                //std::cout<<"I am an A hit first"<<detNum[i]<<std::endl;
                if (detNum[i] == detNum[j]+100){ //this first hit is if the first hit is channel A
                    coinCounter = coinCounter+1;
                }else if (detNum[i]>299){
                    //std::cout<<"I am a B hit first"<<detNum[i]<<std::endl;
                    if (detNum[i] == detNum[j]-100){//this is if the first hit is in B
                        coinCounter = coinCounter+1;
                    }
                }
            }
        }
    }
    //std::cout<<"First Multi= "<<multi<<std::endl;
    for (Int_t j = 0; j < multi; ++j) {
        if (j >= kMaxMulti || detNum[j] < 200 || detNum[j] > 299) continue;
        if (trinCoin == 1){
            if (coinCounter<1) continue; //enforces that a coincidence was seen between A and B channel
        }
        gaggID = detNum[j] - 200;
        if (gaggID < 0 || gaggID >= NGAGG) continue;

        double sumAB = qdc[j][0] + qdc[j][1];
        peak[gaggID]     = qdc[j][3] - (20.0 / 60.0) * sumAB;
        tail[gaggID]     = qdc[j][5] - (55.0 / 60.0) * sumAB;
        traceSum[gaggID] = qdc[j][2] + qdc[j][3] + qdc[j][4] + qdc[j][5] + qdc[j][6] - (115.0 / 60.0) * sumAB;

        if(gaggID>=0 && gaggID<14){//This is here trying to make the ring 1 and 2 PID better. In this loop I will shift the qdc over one.
            peak[gaggID]     = qdc[j][4] - (20.0 / 60.0) * sumAB;
            tail[gaggID]     = qdc[j][6] - (55.0 / 60.0) * sumAB;
            traceSum[gaggID] = qdc[j][3] + qdc[j][4] + qdc[j][5] + qdc[j][6] + qdc[j][7] - (115.0 / 60.0) * sumAB;
        }
        double pidX = 4000.0 * tail[gaggID] / peak[gaggID];
        double pidY = traceSum[gaggID];

        if (std::isfinite(pidX) && std::isfinite(pidY) && peak[gaggID] != 0.0) {
            if (gaggID<=49){
                hPID[gaggID]->Fill(pidX, pidY);
            }else if (gaggID>49){
                hPID[gaggID-50]->Fill(pidX,pidY);
            }
        }
    }
}


//Function to fill the raw germanium hists "We likes it RAW and WRIGGLING"
void makeWrigglingGerm(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, Double_t* e, TH1F** RawGeEnergy){
    for (Int_t j = 0; j < multi; ++j) {
        //std::cout<<"multi"<<multi<<std::endl;
        if (j >= kMaxMulti || detNum[j] < 0 || detNum[j] >= nCrystals) continue;

        RawGeEnergy[detNum[j]]->Fill(e[j]);
    }

}
//Function to fill the BGOs
void makeBGO(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, Double_t* e, TH1F** BGOEnergy){
    for (Int_t j = 0; j < multi; ++j) {
        //std::cout<<"multi"<<multi<<std::endl;
        if (j >= kMaxMulti || (detNum[j]-100) < 0 || (detNum[j]-100) >= nClovers) continue;

        BGOEnergy[(detNum[j]-100)]->Fill(e[j]);
    }

}
//Function to fill the Zero Degrees
void makeZeroDeg(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, Double_t* e, TH1F** ZeroDegEnergy){
    // Detector setup
    const Int_t nZeroDeg  = 2;

    for (Int_t j = 0; j < multi; ++j) {
        //std::cout<<"multi"<<multi<<std::endl;
        if (j >= kMaxMulti || (detNum[j]-400) < 0 || (detNum[j]-400) >= nZeroDeg) continue;

        ZeroDegEnergy[(detNum[j]-400)]->Fill(e[j]);
    }

}
//This makes a a compton suppressed gamma spectrum per crystal
void makeCompSupCalGeEnergy(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, ULong64_t* e_t, Double_t* e, TH1F** CompSupCalGeEnergy, Int_t BGOGeTimeDiff_ns){
    //loop over each hit in the event to check if any are bgo hits
    Double_t calEnergy2 = 0;
    Long64_t timeDiff6 = 0;
    for (Int_t j = 0; j < multi; ++j){
        if(j<kMaxMulti && detNum[j]>99 && detNum[j]<200){
            bgoPresent = true;
        }else{
            bgoPresent = false;
        }   
    }
    //After that loop I know if there is a bgo hit in this event
    for (Int_t j = 0; j < multi; ++j) {
        //std::cout<<"multi"<<multi<<std::endl;
        if (j >= kMaxMulti || detNum[j] < 0 || detNum[j] >= nCrystals) continue; //get Ge hit
        if (bgoPresent){//if there is a BGO hit in this event, loop through and find deltaTs
            for (Int_t k = 0; k < multi; ++k){
                if (k >= kMaxMulti || detNum[k] < 99 || detNum[k] >= 200) continue; //make sure k is a BGO hit'
                timeDiff6 = 10*(static_cast<Long64_t>(e_t[k]) - static_cast<Long64_t>(e_t[j]));
                if(std::abs(timeDiff6)<BGOGeTimeDiff_ns) continue;//This if is actally checking the time. If the BGO and Ge hits are outside of the given window, reject them??? I need to double check logic here.
                calEnergy2 = ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]];
                CompSupCalGeEnergy[detNum[j]]->Fill(calEnergy2);
                }
            }
        if (!bgoPresent){
            calEnergy2 = ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]];
            CompSupCalGeEnergy[detNum[j]]->Fill(calEnergy2);
        }
    }
}

//This makes a compton suppress gamma gamma matrix across all clovers
void makeCompSupCalGammaGamma(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, ULong64_t* e_t, Double_t* e, TH2F* CompSupGammaGamma, Int_t* cloverNum, Int_t BGOGeTimeDiff_ns){
    Double_t x = 0;
    Double_t y = 0;
    bool* check = new bool[kMaxMulti]();
    bgoPresent = false;

    //check if event has bgo hit
    for (Int_t j = 0; j < multi; ++j){
        //std::cout<<"I am here***************************** 1111111111111111111111   "<<detNum[j]<<std::endl;
        if(j<kMaxMulti && detNum[j]>99 && detNum[j]<200){
            bgoPresent = true;
            break;
        }
    }

    if(bgoPresent){//These first nested loops are fill check. Check is true if the timeDiff for that gamma and the BGO is out of the prompt window
        for (Int_t j = 0; j < multi; ++j) {
           //std::cout<<"I am here***************************** 22222222222222222   "<<detNum[j]<<std::endl;
            if (j >= kMaxMulti || detNum[j] < 0 || detNum[j] >= nCrystals) continue;//check that I have a valid Ge hit
            for (Int_t i = 0; i <multi; i++){//loop again, looking for bgo hit
                if(i == j) continue;
                if(detNum[i]>99 && detNum[i]<200){//i is a BGO hit
                    timeDiff_glob = 10*(static_cast<Long64_t>(e_t[i]) - static_cast<Long64_t>(e_t[j]));
                    if(std::abs(timeDiff_glob)>BGOGeTimeDiff_ns){//checking to see if the time diff is outside of the prompt window
                        check[j] = true;
                    }
                }
            }
        }
    

        for (Int_t j = 0; j < multi; ++j) {
            if (j >= kMaxMulti || detNum[j] < 0 || detNum[j] >= nCrystals) continue;
            if(check[j] == false) continue; // skip it if the time diff is in the prompt
            if(e[j] <= 50) continue; //low energy threshold applied
            //std::cout<<"I am here***************************** 333333333333   "<<detNum[j]<<std::endl;
            x = ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]];
            for (Int_t k = 0; k < multi; k++){
                if (k==j || k >= kMaxMulti || detNum[k] < 0 || detNum[k] >= nCrystals) continue; //check for a valid Ge hit AND that it is a different hit than j
                if(check[k] == false) continue;
                if (e[k]<=50) continue;
                y = ge_slope[detNum[k]] * e[k] + ge_intercept[detNum[k]];
                CompSupGammaGamma->Fill(x,y);
            }
        }
    }
    
    if(!bgoPresent){ //if there is not a bgo hit here, we don't worry about vetoing stuff
        for (Int_t j = 0; j < multi; ++j) {
            if (j >= kMaxMulti || detNum[j] < 0 || detNum[j] >= nCrystals) continue;
            if(e[j] <= 50) continue; //low energy threshold applied
            x = ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]];
            for (Int_t k = 0; k < multi; k++){
                if (k==j || k >= kMaxMulti || detNum[k] < 0 || detNum[k] >= nCrystals) continue; //check for a valid Ge hit AND that it is a different hit than j
                if (e[k]<=50) continue;
                y = ge_slope[detNum[k]] * e[k] + ge_intercept[detNum[k]];
                CompSupGammaGamma->Fill(x,y);
            }
        }
    }
    delete[] check;
}
//This will make a deltaT gated gamma gamma
void makeCalGammaGamma(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, ULong64_t* e_t, Double_t* e, TH2F* TotalGammaGamma, Int_t BGOGeTimeDiff_ns){
    Double_t x2 = 0;
    Double_t y2 = 0;
    

    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit
            x2 = ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]];
            for (Int_t k = 0; k < multi; ++k) {//loop over all hits again
                if (detNum[k]>=0 && detNum[k]<nCrystals){//check that a different hit (k=!j) is a Ge hit
                    if(k == j) continue;
                    y2 = ge_slope[detNum[k]] * e[k] + ge_intercept[detNum[k]];
                    if(e[j]>50 && e[k]>50){//Enforcing an energy threshold. A lot of the thresholds were too low and messing with my timing plot
                        TotalGammaGamma->Fill(x2,y2);
                    }
                }
            }
        }
    }
}

void makeAddBackGammaGamma(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, ULong64_t* e_t, Double_t* e, TH2F* AddbackGammaGamma, Int_t BGOGETimeDiff_ns, Int_t* cloverNum, Bool_t* pileup){
        /*loop over events looking for 2 gammas that happen within a prompt window (100 ns)
    Add their energy as x, store their hit number (j and k), then check if multi>2.
    If mult>2, loop through the rest of the events and plot the add back as x and that energy as y. But what if all three are in the time window????
    */

    //***************I think I should fill mult 2 events!!!************************* */
    
    Double_t x3 = 0;
    Double_t y3 = 0;
    Double_t addBackE2 = 0;
    bool* BGOcheck = new bool[kMaxMulti]();
    bool* addBackCheck = new bool[kMaxMulti]();
    bool validAddback = false;
    pairNum = 0;
    Int_t germCount = 0;
    

    timeDiff_glob = 0;

    //Reject pileups
    for (Int_t j = 0; j < multi; ++j) {
        if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit
            germCount = germCount + 1;
        }
        if(pileup[j]){
            return;
        }
    }

    //checking for compton suppression
    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        validAddback = false;
        if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit
            if(e[j]<75) continue;
            for(Int_t k = 0; k<multi; k++){
                if(k == j) continue;
                if (k < kMaxMulti && detNum[k]>99 && detNum[k]<200){//get a BGO hit
                    timeDiff_glob = 10*(static_cast<Long64_t>(e_t[k]) - static_cast<Long64_t>(e_t[j]));
                    if(std::abs(timeDiff_glob)<BGOGETimeDiff_ns){
                        BGOcheck[j] = true;
                    }
                }
            }
        }
    }


    
    if(multi >= 2){
        for (Int_t j = 0; j < multi; ++j) {//loop over all hits
            if(addBackCheck[j]) continue;//this checks to see if this hit has already been used in a addback.
            if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit
                if(e[j]<250) continue; //Low energy threshold, depending on calibration this is (if set to 250) 30 keV or 75 keV
                if(BGOcheck[j] == true) continue;//this gamma was part of a compton scatter into a bgo skip it
                for (Int_t k = 0; k < multi; k++){
                    if (k < kMaxMulti && detNum[k]>=0 && detNum[k]<nCrystals){//get a Ge hit
                        if(k == j) continue;
                        if(addBackCheck[k]) continue;
                        if(e[k]<250) continue;
                        timeDiff_glob = 10*(static_cast<Long64_t>(e_t[k]) - static_cast<Long64_t>(e_t[j]));
                        if(cloverNum[detNum[j]] == cloverNum[detNum[k]]){//here we want the prompt gammas and gammas from teh same clover. Not dealing with diagonals rn
                            pairNum = makeGeTimingPairs(geTimingPairFile, detNum[j], detNum[k]);
                            if(std::abs(timeDiff_glob) <= geTimeGates[pairNum]){
                                addBackE2 = (ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]]) + (ge_slope[detNum[k]] * e[k] + ge_intercept[detNum[k]]);
                                addBackCheck[j] = true;
                                addBackCheck[k] = true;
                                validAddback = true;
                            }
                        }// I still want to keep any gammas that are not addback but are valid gamma gammas
                        if(germCount >= 2 && validAddback == false){
                            x3 = ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]];
                            y3 = ge_slope[detNum[k]] * e[k] + ge_intercept[detNum[k]];
                            AddbackGammaGamma->Fill(x3,y3);

                        }
                        if(multi>2 && validAddback == true){
                            for(Int_t i = 0; i<multi; i++){
                                if(i == j || i == k) continue;
                                if(addBackCheck[i]) continue;
                                if(e[i]<250) continue;
                                if(i < kMaxMulti && detNum[i]>=0 && detNum[i]<nCrystals){//get a Ge hit
                                    y3 = ge_slope[detNum[i]] * e[i] + ge_intercept[detNum[i]];
                                    if(y3 <= 0 || addBackE2 <= 0) continue;
                                    AddbackGammaGamma->Fill(addBackE2,y3);
                                    AddbackGammaGamma->Fill(y3,addBackE2); //making it symmetric
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    delete[] BGOcheck;
    delete[] addBackCheck;
}

ULong64_t makeTrueTimeStamp(Int_t detNum, Int_t cfd, ULong64_t e_t){ // Here I am only passing in the current values of these numbers. This function will NOT loop through the event. Only calculate the good timestamp for the current hit!
    bool boardType100MHz = false;
    bool forcedTriggerBit100 = false;
    ULong64_t forcedTriggerBit500 = 0;
    ULong64_t cfdFractionalTime = 0;
    ULong64_t trueTimeStamp = 0;
    //ULong64_t dumbBigTime = e_t*(4.294967296e9); // this is the timestamp but multiplied by 2^32
    
    if(detNum>=0 && detNum<500){
        boardType100MHz = true;
    }
    if(detNum>=500 && detNum<600){
        boardType100MHz = false;
    }

    //Now we need to split the cfd value into trigger bit and timestamp. Datablock and evtReader keep the 16 bits of cfd information together, I will split the trigger bit and the time here.
    //Let's handle the split for the 100MHz boards firts. I am using the table 4-4 in the Pixie16 manual.
    if(boardType100MHz){
        forcedTriggerBit100 = (cfd >> 15) & 1; //here I am pulling just the last bit in the 16 bit cfd variable.
        cfdFractionalTime = (cfd) & 0x3FFF; // pulling the first 14 bits here. 0x3FFF means 0011 1111 1111 1111.
        if(!forcedTriggerBit100){
            trueTimeStamp = (e_t + (cfdFractionalTime / 32768)) * 10; // calculating true time stamp from the manual page 74
        }
        if(forcedTriggerBit100){
            trueTimeStamp = e_t*10; // there is a 10 ns unit on all of this
        }
    }
    //Now the 500MHz board found in table 4-6
    if(!boardType100MHz){
        forcedTriggerBit500 = (cfd >> 13) & 0x7;
        cfdFractionalTime = (cfd) & 0x1FFF;
        if(forcedTriggerBit500 == 7){//According to the manual, a forced trigger bit is when these three bits are 111. This number is equal to 7 in binary
            trueTimeStamp = 123456789.0;//(e_t*10);
            //std::cout<<"FORCED TRIGGER BIT!!!!"<<std::endl;
        }else{
            trueTimeStamp = ((e_t*50.0) + (forcedTriggerBit500*10.0) - 10 + (cfdFractionalTime / 8192.0)*10.0) * 2.0; // calcualting time stamp with eq 4-7 on page 75
            //std::cout<<"cfd source bits "<<forcedTriggerBit500<<std::endl;
        }
    }
    return trueTimeStamp;
}

void makeCFDSourceBits(Int_t multi, Int_t* detNum, TH1F* forcedTrigger500MHzBoard, Int_t* cfd){
    ULong64_t forcedTriggerBit500 = 0;
    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        if(detNum[j]>=500 && detNum[j]<600){//check for hits in the 500MHz board
            forcedTriggerBit500 = (cfd[j] >> 13) & 0x7;
            forcedTrigger500MHzBoard->Fill(forcedTriggerBit500);
        }
    }
}

void makeGeTimeDiff(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, ULong64_t* e_t, Double_t* e, TH1F** GeTimeDiff, Int_t* cloverNum, Bool_t* pileup){
    Long64_t timeDiff = 0;
    pairNum = 0;
    

    for (Int_t j = 0; j < multi; ++j) {
        if(pileup[j]) return;
    }

    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        //if (pileup[j] == true) return;//skipping the event if it has a pileup flag
        if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit 
            if((ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]])<40) continue;
            for (Int_t k = 0; k < multi; ++k) {//loop over all hits again
                if (k!=j && detNum[k]>=0 && detNum[k]<nCrystals){//check that a different hit (k=!j) is a Ge hit
                    //if (pileup[k] == true) return;
                    if (cloverNum[detNum[j]] == cloverNum[detNum[k]]){//Check that the hits are in the same clover
                        pairNum = makeGeTimingPairs(geTimingPairFile, detNum[j], detNum[k]);
                        //std::cout<<"I made it through the same clover check pairID: "<<pairNum<<std::endl;
                        timeDiff = 10 * (static_cast<Long64_t>(e_t[k]) - static_cast<Long64_t>(e_t[j]));//find the time diff between them
                        if((ge_slope[detNum[k]] * e[k] + ge_intercept[detNum[k]])>40){
                            GeTimeDiff[pairNum]->Fill(timeDiff);// fill that diff the diff hist
                        }
                    }
                }
            }
        }
    }
}
void makeTotalGeTimeDiff(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, ULong64_t* e_t, Double_t* e, TH1F* TotalGeTimeDiff, Int_t* cloverNum){
    Long64_t timeDiff2 = 0;

    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit
            for (Int_t k = 0; k < multi; ++k) {//loop over all hits again
                if (k!=j && detNum[k]>=0 && detNum[k]<nCrystals){//check that a different hit (k=!j) is a Ge hit
                    timeDiff2 = 10*(static_cast<Long64_t>(e_t[k]) - static_cast<Long64_t>(e_t[j]));//find the time diff between them
                    //std::cout<<"j timestamp: "<<detNum[j]<<" | "<<e_t[j]<<" | k timestamp: "<<detNum[k]<<" | "<<e_t[k]<<" | timediff: "<<timeDiff<<std::endl;
                    if(e[j]>50 && e[k]>50){//Enforcing an energy threshold. A lot of the thresholds were too low and messing with my timing plots
                        TotalGeTimeDiff->Fill(timeDiff2);// fill that diff the diff hist
                    }
                }
            }
        }
    }
}
void makeLabrTimeDiff(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, TH1F** LabrTimeDiff, Bool_t* pileup, Int_t* cfd, UInt_t (*qdc)[8], ULong64_t* e_t){
    Long64_t timeDiff = 0;
    pairNum = 0;
    Double_t labrQDCEnergy1;
    Double_t labrQDCEnergy2;
    ULong64_t time1;
    ULong64_t time2;
    

    //for (Int_t j = 0; j < multi; ++j) {
        //if(pileup[j]) return;      
    //}


    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        labrQDCEnergy1 = ((qdc[j][2]+qdc[j][3])-qdc[j][1])*(100.0/140.0);
        labrQDCEnergy1 = labrQDCEnergy1*labr_slope[detNum[j]-500]+labr_intercept[detNum[j]-500];

        //if (pileup[j] == true) return;//skipping the event if it has a pileup flag
        if (j < kMaxMulti && detNum[j]>=500 && detNum[j]<600){//get a labr hit 
            if(labrQDCEnergy1>=1000 && labrQDCEnergy1<=2100){
                time1 = makeTrueTimeStamp(detNum[j], cfd[j], e_t[j]);
                if(time1 == 123456789){
                    //std::cout<<"Forced bit in t1*******************"<<std::endl;
                    continue;
                }
                for (Int_t k = 0; k < multi; ++k) {//loop over all hits again
                    if (k!=j && detNum[k]>=500 && detNum[k]<600){//check that a different hit (k=!j) is a hit
                        if (detNum[j] == detNum[k]) continue; // check if those two hits are the same detector firing twice.
                        time2 = makeTrueTimeStamp(detNum[k], cfd[k], e_t[k]);
                        labrQDCEnergy2 = ((qdc[k][2]+qdc[k][3]) - qdc[k][1])*(100.0/140.0);
                        labrQDCEnergy2 = labrQDCEnergy2*labr_slope[detNum[k]-500]+labr_intercept[detNum[k]-500];
                        if(time2 == 123456789){
                            //std::cout<<"Forced bit in t2************"<<std::endl;
                            continue;
                        }
                        if((labrQDCEnergy2>=1000 && labrQDCEnergy2<= 2100)){
                            pairNum = makeLabrTimingPairs(labrTimingPairFile, (detNum[j] - 500), (detNum[k] - 500));
                           // std::cout<<"I made it through energy2 check "<<pairNum<<std::endl;
                            timeDiff = time1 - time2;//find the time diff between them
                            LabrTimeDiff[pairNum]->Fill(abs(timeDiff));// fill that diff the diff hist
                            //std::cout<<"I am in labr time diff fill statement************************************"<<std::endl;
                        }
                    }
                }
            }
        }
    }
}

// void makeGAGGGAGGTimeDiff(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, TH1F** GAGGGAGGTimeDiff, Bool_t* pileup, Int_t* cfd, UInt_t (*qdc)[8], ULong64_t* e_t){
//     Long64_t timeDiff = 0;
//     pairNum = 0;
    
// ///THIS NEEDS TO HAPPEN!!!! I HAVE NOT FINISHED THIS SCRIPT!!! I AM PAUSING HERE TO FIGURE OUT TRACES
//     //for (Int_t j = 0; j < multi; ++j) {
//         //if(pileup[j]) return;      
//     //}


//     for (Int_t j = 0; j < multi; ++j) {//loop over all hits
//         //if (pileup[j] == true) return;//skipping the event if it has a pileup flag
//         if (j < kMaxMulti && detNum[j]>=500 && detNum[j]<600){//get a Ge hit 
//             //if((labr_slope[detNum[j] - 450] * labrQDCEnergy1 + labr_intercept[detNum[j] - 450])<1) continue; //the -450 here is for the calibration read mapping. Labr start at index 50 there. Therefore 500 - 450. 
//             for (Int_t k = 0; k < multi; ++k) {//loop over all hits again
//                 if (k!=j && detNum[k]>=500 && detNum[k]<600){//check that a different hit (k=!j) is a Ge hit
//                     if (detNum[j] == detNum[k]) continue; // check if those two hits are the same detector firing twice. 
//                         pairNum = makeLabrTimingPairs(labrTimingPairFile, (detNum[j] - 500), (detNum[k] - 500));
//                         //std::cout<<"I made it through the same clover check pairID: "<<pairNum<<std::endl;
//                         timeDiff = (makeTrueTimeStamp(detNum[k], cfd[k], e_t[k]) - makeTrueTimeStamp(detNum[j], cfd[j], e_t[j]));//find the time diff between them
//                         //if((labr_slope[detNum[k] - 450] * labrQDCEnergy2 + labr_intercept[detNum[k] - 450])>1){
//                             LabrTimeDiff[pairNum]->Fill(abs(timeDiff));// fill that diff the diff hist
//                             //std::cout<<"I am in labr time diff fill statement************************************"<<std::endl;
//                         //}
//                 }
//             }
//         }
//     }
// }

void makeTotalLaBr3TimeDiff(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, ULong64_t* e_t, Double_t* e, TH1F* TotalLaBr3TimeDiff, Int_t* cloverNum, Int_t* cfd){
    Double_t timeDiff3 = 0;

    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        if (j < kMaxMulti && detNum[j]>=500 && detNum[j]<600){//get a Labr hit
            //std::cout<<"I GOT IN THE FIRST LOOP"<<std::endl;
            for (Int_t k = 0; k < multi; ++k) {//loop over all hits again
                if (k!=j && detNum[k]>=500 && detNum[k]<600){//check that a different hit (k=!j) is a LabR3 hit
                    timeDiff3 = makeTrueTimeStamp(detNum[j], cfd[j], e_t[j]) - makeTrueTimeStamp(detNum[k], cfd[k], e_t[k]); // this has already been converted to ns!!
                    if(timeDiff3 == 0 || std::fmod(timeDiff3,10.0) == 0){
                        continue;
                    }
                    //std::cout<<"Calc'd time diff "<<timeDiff3<<std::endl;
                    //timeDiff3 = 10*(static_cast<Long64_t>(e_t[k]) - static_cast<Long64_t>(e_t[j]));//find the time diff between them
                    //std::cout<<"j timestamp: "<<detNum[j]<<" | "<<e_t[j]<<" | k timestamp: "<<detNum[k]<<" | "<<e_t[k]<<" | timediff: "<<timeDiff<<std::endl;
                    //if(e[j]>50 && e[k]>50){//Enforcing an energy threshold. A lot of the thresholds were too low and messing with my timing plots
                        TotalLaBr3TimeDiff->Fill(timeDiff3);// fill that diff the diff hist
                    //}
                }
            }
        }
    }
}

void makeLabrTimeDiffvLabrEnergy(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, TH2F** LabrTimeDiffvLabrEnergy, Bool_t* pileup, Int_t* cfd, UInt_t (*qdc)[8], ULong64_t* e_t){
    Long64_t timeDiff = 0;
    pairNum = 0;
    Double_t labrQDCEnergy1;
    Double_t labrQDCEnergy2;
    ULong64_t time1;
    ULong64_t time2;
    

    //for (Int_t j = 0; j < multi; ++j) {
        //if(pileup[j]) return;      
    //}


    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        labrQDCEnergy1 = ((qdc[j][2]+qdc[j][3])-qdc[j][1])*(100.0/140.0);
        labrQDCEnergy1 = labrQDCEnergy1*labr_slope[detNum[j]-500]+labr_intercept[detNum[j]-500];

        //if (pileup[j] == true) return;//skipping the event if it has a pileup flag
        if (j < kMaxMulti && detNum[j]>=500 && detNum[j]<600){//get a labr hit 
            if(labrQDCEnergy1>=0 && labrQDCEnergy1<=5000){
                time1 = makeTrueTimeStamp(detNum[j], cfd[j], e_t[j]);
                if(time1 == 123456789){
                    //std::cout<<"Forced bit in t1*******************"<<std::endl;
                    continue;
                }
                for (Int_t k = 0; k < multi; ++k) {//loop over all hits again
                    if (k!=j && detNum[k]>=500 && detNum[k]<600){//check that a different hit (k=!j) is a hit
                        if (detNum[j] == detNum[k]) continue; // check if those two hits are the same detector firing twice.
                        time2 = makeTrueTimeStamp(detNum[k], cfd[k], e_t[k]);
                        labrQDCEnergy2 = ((qdc[k][2]+qdc[k][3]) - qdc[k][1])*(100.0/140.0);
                        labrQDCEnergy2 = labrQDCEnergy2*labr_slope[detNum[k]-500]+labr_intercept[detNum[k]-500];
                        if(time2 == 123456789){
                            //std::cout<<"Forced bit in t2************"<<std::endl;
                            continue;
                        }
                        if((labrQDCEnergy2>=0 && labrQDCEnergy2<= 5000)){
                            pairNum = makeLabrTimingPairs(labrTimingPairFile, (detNum[j] - 500), (detNum[k] - 500));
                           // std::cout<<"I made it through energy2 check "<<pairNum<<std::endl;
                            timeDiff = time1 - time2;//find the time diff between them
                            LabrTimeDiffvLabrEnergy[pairNum]->Fill(abs(timeDiff),labrQDCEnergy1);// fill that diff the diff hist
                            //std::cout<<"I am in labr time diff fill statement************************************"<<std::endl;
                        }
                    }
                }
            }
        }
    }
}

void makeBGOGeTimeDiff(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, ULong64_t* e_t, Double_t* e, TH1F** BGOGeTimeDiff, Int_t* cloverNum){
    Long64_t timeDiff5 = 0;

    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit
            for (Int_t k = 0; k < multi; ++k) {//loop over all hits again
                if (k!=j && detNum[k]>=100 && detNum[k]<200){//check that a different hit (k=!j) is a BGO hit
                    timeDiff5 = 10*(static_cast<Long64_t>(e_t[k]) - static_cast<Long64_t>(e_t[j]));//find the time diff between them
                    if(e[j]>200){//Enforcing an energy threshold. A lot of the thresholds were too low and messing with my timing plots
                        if(cloverNum[detNum[j]] == (detNum[k]-100)){
                            BGOGeTimeDiff[detNum[j]]->Fill(timeDiff5);//std::abs(timeDiff5));// fill that diff the diff hist
                        }
                    }
                }
            }
        }
    }
}
void readAlphaCutsFile(const char* alphaCutFile, TCutG** alphaCuts) {
    // Open the ROOT file
    TFile* file = TFile::Open(alphaCutFile, "READ");
    if (!file || file->IsZombie()) {
        std::cerr << "Error opening cut file: " << alphaCutFile << std::endl;
        return;
    }

    // Regular expression to match "hPID##_alpha"
    std::regex pattern("PID(\\d{2})_alpha");

    // Loop through all keys in the file
    TIter next(file->GetListOfKeys());
    TKey* key;
    while ((key = (TKey*)next())) {
        TObject* obj = key->ReadObj();

        // Only proceed if object is a TCutG
        if (!obj->InheritsFrom("TCutG")) continue;

        std::string name = obj->GetName();
        std::smatch match;
        if (std::regex_match(name, match, pattern)) {
            int idx = std::stoi(match[1]);
            if (idx >= 0 && idx < 64) {
                alphaCuts[idx] = (TCutG*)obj;
                std::cout << "Loaded alpha cut for detector " << idx << std::endl;
            } else {
                std::cerr << "Warning: Cut index " << idx << " is out of bounds (0-63)" << std::endl;
            }
        } else {
            std::cerr << "Skipping unrecognized cut: " << name << std::endl;
        }
    }

    file->Close();
}



void readCalFile(const char* calFile){
    std::ifstream infile(calFile);
    if (!infile.is_open()) {
        std::cerr << "Error: Could not open calibration file " << calFile << std::endl;
        return;
    }

    int detID;
    double slope, intercept;

    // Initialize all entries to default values
    for (int i = 0; i < nCrystals + nBEGe + nLaBr3; ++i) {
        ge_slope[i] = 1.0;      // default: slope=1 (no gain correction)
        ge_intercept[i] = 0.0;   // default: intercept=0 (no offset correction)
    }

    // Read the file line by line
    std::string line;
    while (std::getline(infile, line)) {
        std::istringstream iss(line);
        if (!(iss >> detID >> slope >> intercept)) {
            std::cerr << "Warning: Bad line in calibration file: " << line << std::endl;
            continue;
        }
        if (detID >= 0 && detID < nCrystals + nBEGe) {
            ge_slope[detID] = slope;
            ge_intercept[detID] = intercept;
        } 
        if (detID>nCrystals + nBEGe){
            labr_slope[detID - 50] = slope;
            labr_intercept[detID- 50] = intercept;
        } 
        
        if (detID> nCrystals + nBEGe + nLaBr3 + 5) {
            std::cerr << "Warning: Detector ID out of range: " << detID << std::endl;
        }
    }

    infile.close();

    infile.close();
    std::cout << "Calibration loaded successfully from " << calFile << std::endl;
}

void readGeTimeGatesFile(const char* timeGatesFile){
    std::ifstream infile(timeGatesFile);
    if (!infile.is_open()) {
        std::cerr << "Error: Could not open time gates file " << timeGatesFile << std::endl;
        return;
    }

    int detID2;
    double timeGateRead;

    // Initialize all entries to default values
    for (int i = 0; i < nCrystalPairs; ++i) {
        geTimeGates[i] = 1.0;      // default: slope=1 (no gain correction)
    }

    // Read the file line by line
    std::string line;
    while (std::getline(infile, line)) {
        std::istringstream iss(line);
        if (!(iss >> detID2 >> timeGateRead)) {
            std::cerr << "Warning: Bad line in time gates file: " << line << std::endl;
            continue;
        }
        if (detID2 >= 0 && detID2 < nCrystalPairs) {
            geTimeGates[detID2] = timeGateRead;
        } else {
            std::cerr << "Warning: Detector ID out of range: " << detID2 << std::endl;
        }
    }

    infile.close();

    infile.close();
    std::cout << "Prompt time gates loaded successfully from " << timeGatesFile << std::endl;
}
void makeCalGeEnergy(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, Double_t* e, TH1F** CalGeEnergy){
    //Make a varible to hold the calEnergy
    Double_t calEnergy = 0;

    for (Int_t j = 0; j < multi; ++j) {
        //std::cout<<"multi"<<multi<<std::endl;
        if (j >= kMaxMulti || detNum[j] < 0 || detNum[j] >= nCrystals) continue;
        //Apply a linear calibration 
        calEnergy = ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]];
        CalGeEnergy[detNum[j]]->Fill(calEnergy);
    }

}
void makeCloverNum(Int_t* cloverNum){
    for (Int_t j = 0; j<nCrystals; j++){
        cloverNum[j]=j/4;
    }
}



void makeAddBackGeEnergy(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, Double_t* e, Int_t* cloverNum, TH1F** AddBackGeEnergy, bool includeDiagonals){//I am only doing 2 hits here.
    Double_t addBackEnergy = 0;
    Double_t singleCalEnergy = 0;
    //This makes all of the events that I addback are compton suppressed in the sense that there are no BGO hits.
    //for (Int_t j = 0; j < multi; ++j){
        //if(j<kMaxMulti && detNum[j]>99 && detNum[j]<200)
            //return;
    //}
    for(Int_t j = 0; j<multi; j++){
        if (j >= kMaxMulti || detNum[j] < 0 || detNum[j] >= nCrystals) continue; //this was not a Ge hit
        if (multi == 1){
            singleCalEnergy = (ge_slope[detNum[j]]*e[j]) + ge_intercept[detNum[j]];
            AddBackGeEnergy[cloverNum[detNum[j]]]->Fill(singleCalEnergy);
        }else{
            for(Int_t k = j+1; k<multi; k++){//starting this at j+1 to avoid double counting.
                if (k >= kMaxMulti || detNum[k] < 0 || detNum[k] >= nCrystals || k == j) continue; //this was not a Ge hit or it is the same hit

                if (cloverNum[detNum[k]] == cloverNum[detNum[j]]){//Are the hits in the same clover
                    if (includeDiagonals == true){
                        addBackEnergy = (ge_slope[detNum[j]]*e[j]+ge_intercept[detNum[j]])+(ge_slope[detNum[k]]*e[k]+ge_intercept[detNum[k]]);
                        AddBackGeEnergy[cloverNum[detNum[j]]]->Fill(addBackEnergy); 
                    }else{
                        /*
                        std::cout<<"Multiplicity: "<<multi<<" det number one: "<<detNum[j]<<" Hit one energy: "<<e[j]<<std::endl;
                        std::cout<<"Multiplicity: "<<multi<<" Det number two: "<<detNum[k]<<" Hit two energy: "<< e[k] <<std::endl;
                        std::cout<<"Slope one: "<<ge_slope[detNum[j]]<<" Slope two: "<<ge_slope[detNum[k]]<<std::endl;
                        std::cout<<"__________________________________________________________________________________________________"<<std::endl;
                        
                        addBackEnergy = ((ge_slope[detNum[j]]*e[j])+ge_intercept[detNum[j]])+((ge_slope[detNum[k]]*e[k])+ge_intercept[detNum[k]]);
                        AddBackGeEnergy[cloverNum[detNum[j]]]->Fill(addBackEnergy);
                        */
                        
                        //Below is logic for excluding the diagonals.
                        
                        if (detNum[j]%4 == 0 || detNum[j]%4 == 2){//Is the first hit a blue or green hit?
                            if(detNum[k]%4 == 1 || detNum[k]%4 == 3){//Is the second hit a red or black hit? I only want to add hits that happened in adjacent crystals. No diagonals
                                //std::cout<<"Scenario A | "<<"Multiplicity: "<<multi<< " det number one: "<<detNum[j]<< " Hit one energy: "<<e[j]<<" det number two: "<<detNum[k]<<" Hit two energy: "<< e[k] <<" Clover number one: "<<cloverNum[detNum[j]]<<" Clover number two: "<<cloverNum[detNum[k]]<<std::endl;  
                                addBackEnergy = (ge_slope[detNum[j]]*e[j]+ge_intercept[detNum[j]])+(ge_slope[detNum[k]]*e[k]+ge_intercept[detNum[k]]);
                                AddBackGeEnergy[cloverNum[detNum[j]]]->Fill(addBackEnergy); 
                            }
                        }
                        if (detNum[j]%4 == 1 || detNum[j]%4 == 3){//Vice versa of above
                            if(detNum[k]%4 == 0 || detNum[k]%4 == 2){
                                //std::cout<<"Scenario B | "<<"Multiplicity: "<<multi<< " det number one: "<<detNum[j]<< " Hit one energy: "<<e[j]<<" det number two: "<<detNum[k]<<" Hit two energy: "<< e[k] <<" Clover number one: "<<cloverNum[detNum[j]]<<" Clover number two: "<<cloverNum[detNum[k]]<<std::endl;
                                addBackEnergy = (ge_slope[detNum[j]]*e[j]+ge_intercept[detNum[j]])+(ge_slope[detNum[k]]*e[k]+ge_intercept[detNum[k]]);
                                AddBackGeEnergy[cloverNum[detNum[j]]]->Fill(addBackEnergy); 
                            }
                        }
                    }
                }
            }
        }
    }
}

void makeCompSupCloverSumEnergy(TH1F** CompSupCalGeEnergy, TH1F** CompSupCloverSumEnergy){
    for (int i = 0; i < 11; ++i) {
        int idx0 = i * 4;
        CompSupCloverSumEnergy[i]->Reset(); // clear existing bin contents

        for (int j = 0; j < 4; ++j) {
            CompSupCloverSumEnergy[i]->Add(CompSupCalGeEnergy[idx0 + j]);
        }
    }

}

void makeAlphaCutPID(TH2F** hPID, TCutG** alphaCuts, TH2F** hPID_alpha) {
    for (int i = 0; i < NGAGG; ++i) {
        if (!hPID[i] || !alphaCuts[i] || !hPID_alpha[i]) continue;  // Skip if any pointer is null
        //std::cout<<"I GOT HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<std::endl;
        hPID_alpha[i]->Reset();  // Ensure we start with an empty histogram

        int nBinsX = hPID[i]->GetNbinsX();
        int nBinsY = hPID[i]->GetNbinsY();

        for (int xBin = 1; xBin <= nBinsX; ++xBin) {
            double x = hPID[i]->GetXaxis()->GetBinCenter(xBin);
            for (int yBin = 1; yBin <= nBinsY; ++yBin) {
                double y = hPID[i]->GetYaxis()->GetBinCenter(yBin);
                if (alphaCuts[i]->IsInside(x, y)) {
                    double content = hPID[i]->GetBinContent(xBin, yBin);
                    hPID_alpha[i]->SetBinContent(xBin, yBin, content);
                    //std::cout<<"I GOT HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<std::endl;
                }
            }
        }
    }
}

bool checkAlphaPID(TTree* tree, int kMaxMulti, int NGAGG, Int_t multi, Int_t* detNum, UInt_t (*qdc)[8], TCutG** alphaCuts){
    //int trinMult = 0;
    int numAlphaHits = 0;
    for (Int_t j = 0; j < multi; ++j) {
        if (j >= kMaxMulti || detNum[j] < 200 || detNum[j] > 299) continue; //Usiung this to get a Trinity hit
        //Calc pidX and Y
        gaggID = detNum[j] - 200;
        if (gaggID < 0 || gaggID >= NGAGG) continue;

        double sumAB2 = qdc[j][0] + qdc[j][1];
        peak[gaggID]     = qdc[j][3] - (20.0 / 60.0) * sumAB2;
        tail[gaggID]     = qdc[j][5] - (55.0 / 60.0) * sumAB2;
        traceSum[gaggID] = qdc[j][2] + qdc[j][3] + qdc[j][4] + qdc[j][5] + qdc[j][6] - (115.0 / 60.0) * sumAB2;

        if(gaggID>=0 && gaggID<14){//This is here trying to make the ring 1 and 2 PID better. In this loop I will shift the qdc over one.
            peak[gaggID]     = qdc[j][4] - (20.0 / 60.0) * sumAB2;
            tail[gaggID]     = qdc[j][6] - (55.0 / 60.0) * sumAB2;
            traceSum[gaggID] = qdc[j][3] + qdc[j][4] + qdc[j][5] + qdc[j][6] + qdc[j][7] - (115.0 / 60.0) * sumAB2;
        }
        double pidX2 = 4000.0 * tail[gaggID] / peak[gaggID];
        double pidY2 = traceSum[gaggID];

        //Check if it is in the cut
        if(alphaCuts[gaggID] && alphaCuts[gaggID]->IsInside(pidX2,pidY2)){
            numAlphaHits = numAlphaHits + 1;
        }
    }
    //return true if it is, false if not
    if(numAlphaHits >= 1){
        return true;
    }else{
        return false;
    }
    //std::cout<<"Trin Mult: "<<trinMult<<std::endl;
}

void makeAlphaCutAddBackGammaGamma(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, UInt_t (*qdc)[8], ULong64_t* e_t, Double_t* e, TH2F* AlphaAddbackGammaGamma, Int_t BGOGETimeDiff_ns, Int_t* cloverNum, Bool_t* pileup, TCutG** alphaCuts){

    
    Double_t x4 = 0;
    Double_t y4 = 0;
    Double_t addBackE3 = 0;
    bool* BGOcheck2 = new bool[kMaxMulti]();
    bool* addBackCheck2 = new bool[kMaxMulti]();
    bool validAddback2 = false;
    pairNum = 0;
    Int_t germCount2 = 0;
    timeDiff_glob = 0;

    //Check to see if the event has an alpha in it
    alphaCheck = checkAlphaPID(tree, kMaxMulti, NGAGG, multi, detNum, qdc, alphaCuts);

    if(!alphaCheck) return;


    //Reject pileups
    for (Int_t j = 0; j < multi; ++j) {
        if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit
            germCount2 = germCount2 + 1;
        }
        if(pileup[j]){
            return;
        }
    }

    //checking for compton suppression
    for (Int_t j = 0; j < multi; ++j) {//loop over all hits
        validAddback2 = false;
        if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit
            if(e[j]<75) continue;
            for(Int_t k = 0; k<multi; k++){
                if(k == j) continue;
                if (k < kMaxMulti && detNum[k]>99 && detNum[k]<200){//get a BGO hit
                    timeDiff_glob = 10*(static_cast<Long64_t>(e_t[k]) - static_cast<Long64_t>(e_t[j]));
                    if(std::abs(timeDiff_glob)<BGOGETimeDiff_ns){
                        BGOcheck2[j] = true;
                    }
                }
            }
        }
    }


    
    if(multi >= 2){
        for (Int_t j = 0; j < multi; ++j) {//loop over all hits
            if(addBackCheck2[j]) continue;//this checks to see if this hit has already been used in a addback.
            if (j < kMaxMulti && detNum[j]>=0 && detNum[j]<nCrystals){//get a Ge hit
                if(e[j]<250) continue; //Low energy threshold, depending on calibration this is (if set to 250) 30 keV or 75 keV
                if(BGOcheck2[j] == true) continue;//this gamma was part of a compton scatter into a bgo skip it
                for (Int_t k = 0; k < multi; k++){
                    if (k < kMaxMulti && detNum[k]>=0 && detNum[k]<nCrystals){//get a Ge hit
                        if(k == j) continue;
                        if(addBackCheck2[k]) continue;
                        if(e[k]<250) continue;
                        timeDiff_glob = 10*(static_cast<Long64_t>(e_t[k]) - static_cast<Long64_t>(e_t[j]));
                        if(cloverNum[detNum[j]] == cloverNum[detNum[k]]){//here we want the prompt gammas and gammas from teh same clover. Not dealing with diagonals rn
                            pairNum = makeGeTimingPairs(geTimingPairFile, detNum[j], detNum[k]);
                            if(std::abs(timeDiff_glob) <= geTimeGates[pairNum]){
                                addBackE3 = (ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]]) + (ge_slope[detNum[k]] * e[k] + ge_intercept[detNum[k]]);
                                addBackCheck2[j] = true;
                                addBackCheck2[k] = true;
                                validAddback2 = true;
                            }
                        }// I still want to keep any gammas that are not addback but are valid gamma gammas
                        if(germCount2 >= 2 && validAddback2 == false){
                            x4 = ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]];
                            y4 = ge_slope[detNum[k]] * e[k] + ge_intercept[detNum[k]];
                            AlphaAddbackGammaGamma->Fill(x4,y4);

                        }
                        if(multi>2 && validAddback2 == true){
                            for(Int_t i = 0; i<multi; i++){
                                if(i == j || i == k) continue;
                                if(addBackCheck2[i]) continue;
                                if(e[i]<250) continue;
                                if(i < kMaxMulti && detNum[i]>=0 && detNum[i]<nCrystals){//get a Ge hit
                                    y4 = ge_slope[detNum[i]] * e[i] + ge_intercept[detNum[i]];
                                    if(y4 <= 0 || addBackE3 <= 0) continue;
                                    AlphaAddbackGammaGamma->Fill(addBackE3,y4);
                                    AlphaAddbackGammaGamma->Fill(y4,addBackE3); //making it symmetric
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    delete[] BGOcheck2;
    delete[] addBackCheck2;
}

void makeTraces(Int_t multi, Int_t* detNum, Int_t (*trace)[280], TH2F** Traces){
    for (Int_t j = 0; j < multi; ++j){
        if(detNum[j]>=500 && detNum[j]<600){//hard coding in to only look at labr traces
            for(Int_t k = 0; k<280; k++){//looping over trace length hardcoded length for now
                Traces[detNum[j]-500]->Fill(k, trace[j][k]);
                //std::cout<<"j: "<<j<<" Trace[j][k] = "<<trace[j][k]<<std::endl;
            }
        }
    }
}

// Create a TChain from all .root files in a directory
TChain* makeTChain(const char* inputDir, const char* treeName = "tree") {
    TChain* chain = new TChain(treeName);

    void* dirp = gSystem->OpenDirectory(inputDir);
    if (!dirp) {
        std::cerr << "Error: Cannot open directory: " << inputDir << std::endl;
        return nullptr;
    }

    const char* file;
    while ((file = gSystem->GetDirEntry(dirp))) {
        TString fname = file;
        if (!fname.EndsWith(".root")) continue;

        TString fullPath = TString(inputDir) + fname;
        chain->Add(fullPath);
    }

    if (chain->GetEntries() == 0) {
        std::cerr << "Warning: No entries found in TChain from directory: " << inputDir << std::endl;
        delete chain;
        return nullptr;
    }

    return chain;
}

TTree* makeOpenRootTree(const char* filename) {
    TFile* file = TFile::Open(filename, "READ");
    if (!file || file->IsZombie()) {
        std::cerr << "Error opening file: " << filename << std::endl;
        return nullptr;
    }

    TTree* tree = nullptr;
    file->GetObject("tree", tree);   // MUST match tree name

    if (!tree) {
        std::cerr << "Error: TTree 'tree' not found in file: "
                  << filename << std::endl;
        file->ls();
        return nullptr;
    }

    // DO NOT call tree->SetDirectory(nullptr)
    // DO NOT close or delete file

    return tree;
}



//Function to fill the LaBr3s
void makeLaBr3(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, Double_t* e, TH1F** LaBr3Energy, UInt_t (*qdc)[8]){
    for (Int_t j = 0; j < multi; ++j) {
        //std::cout<<"multi"<<multi<<std::endl;

        
        if (j >= kMaxMulti || (detNum[j] - 500) < 0 || (detNum[j] - 500) >= nLaBr3) continue;
        double labr3EnergySum = (((double)qdc[j][2] + (double)qdc[j][3]) - (double)qdc[j][1]) * (100.0/140.0);
        
        labr3EnergySum = labr_slope[detNum[j] - 500] * labr3EnergySum + labr_intercept[detNum[j] - 500];
        //std::cout<<"DetNum - 500: "<<detNum[j] - 500<<" Slope: "<<labr_slope[detNum[j] - 500]<<" Intercept: "<<labr_intercept[detNum[j] - 500]<<std::endl;
        //std::cout<<" Sum energy: "<<labr3EnergySum<<" QDC1: "<<qdc[j][1]<<std::endl;
        LaBr3Energy[(detNum[j]-500)]->Fill(labr3EnergySum);
    }
}

//Function to fill the BEGe channels
void makeRawBEGe(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, Double_t* e, TH1F** RawBEGe){
    for (Int_t j = 0; j < multi; ++j) {
        //std::cout<<"multi"<<multi<<std::endl;
        if (j >= kMaxMulti || (detNum[j]) < 44 || (detNum[j]) > 45) continue;
        //std::cout<<"detNum "<<detNum[j]<<std::endl;
        RawBEGe[(detNum[j]) - 44]->Fill(e[j]);
    }
}

void makeCalBEGe(TTree* tree, int kMaxMulti, Int_t multi, Int_t* detNum, Double_t* e, TH1F** calBEGe){
    //Make a varible to hold the calEnergy
    Double_t calEnergy3 = 0;

    for (Int_t j = 0; j < multi; ++j) {
        //std::cout<<"multi"<<multi<<std::endl;
        if (j >= kMaxMulti || detNum[j] < 44 || detNum[j] > 45) continue;
        //Apply a linear calibration 
        calEnergy3 = ge_slope[detNum[j]] * e[j] + ge_intercept[detNum[j]];
        calBEGe[detNum[j] - 44]->Fill(calEnergy3);
    }
}

void Process_Events(const char* inputFile, 
                    const char* calFilename,
                    const char* gePromptGatesFile,
                    const char* outputDirectory,
                    const char* outputFilename, 
                    bool trinCoin,
                    bool includeDiagonals,
                    const char* alphaCutsFile,
                    const char* geTimeingPairFileRead,
                    Int_t BGOGeTimeGate_ns,
                    const char* labrTimeingPairFileRead) {
    //Read in the calobration File
    readCalFile(calFilename);
    readGeTimeGatesFile(gePromptGatesFile);
    geTimingPairFile = geTimeingPairFileRead;
    labrTimingPairFile = labrTimeingPairFileRead;

    // Open the input ROOT file
    //TChain* tree = makeTChain(inputDir, "tree");  // this chains together all of the root files in a directory
    TTree* tree = makeOpenRootTree(inputFile);
    if(!tree){
        std::cout<<"No Tree Found"<<std::endl;
        return;
    }

    // Set up variables to hold the branch data
    const Int_t kMaxMulti = 500;
    Int_t multi;
    Int_t detNum[kMaxMulti];
    Double_t e[kMaxMulti];
    UInt_t qdc[kMaxMulti][8];
    ULong64_t e_t[kMaxMulti];
    Int_t cloverNum[nCrystals];
    Bool_t pileup[kMaxMulti];
    Int_t cfd[kMaxMulti];
    Int_t trace[416][280]; //hard coded in tracelength = 280 here.

    tree->SetBranchAddress("multi", &multi);
    tree->SetBranchAddress("detID", detNum);
    tree->SetBranchAddress("e", e);
    tree->SetBranchAddress("e_t", e_t);
    tree->SetBranchAddress("pileup", pileup);
    tree->SetBranchAddress("qdc", qdc);
    tree->SetBranchAddress("cfd",cfd);
    tree->SetBranchAddress("trace",trace);

   


    //Make the cloverNum array
    makeCloverNum(cloverNum);
    // Initialize all elements with NaN
    for (int i = 0; i < NGAGG; ++i){
            traceSum[i] = peak[i] = tail[i] = TMath::QuietNaN();
    }
    // Create gamma energy histograms
    TH1F* RawGeEnergy[nCrystals];
    for (int i = 0; i < nCrystals; ++i) {
        RawGeEnergy[i] = new TH1F(Form("RawGeEnergy%02d", i), Form("e - %02d", i), 4000, 0, 16000);
        RawGeEnergy[i]->SetDirectory(nullptr);//comment this out if you want to use this histogram
    }
    //Create bgo energy hists
    TH1F* BGOEnergy[nClovers];
    for (int i = 0; i < nClovers; ++i) {
        BGOEnergy[i] = new TH1F(Form("BGOEnergy%02d", i), Form("e - %02d", i), 20000, 0, 20000);
        BGOEnergy[i]->SetDirectory(nullptr);//comment this out if you want to use this histogram
    }
    //Create zerodeg hists
    TH1F* ZeroDegEnergy[nZeroDeg];
    for (int i = 0; i < nZeroDeg; ++i) {
        ZeroDegEnergy[i] = new TH1F(Form("ZeroDegEnergy%02d", i), Form("e - %02d", i), 20000, 0, 20000);
        ZeroDegEnergy[i]->SetDirectory(nullptr);//comment this out if you want to use this histogram
    }
    //Create the compton suppressed Ge Hists
    TH1F* CompSupCalGeEnergy[nCrystals];
    for (int i = 0; i < nCrystals; ++i) {
        CompSupCalGeEnergy[i] = new TH1F(Form("CompSupCalGeEnergy%02d", i), Form("e - %02d", i), 4000, 0, 4000);
        CompSupCalGeEnergy[i]->SetDirectory(nullptr);//comment this out if you want to use this histogram
    }
    // Create calibrated gamma energy histograms
    TH1F* CalGeEnergy[nCrystals];
    TH1F* BGOGeTimeDiff[nCrystals];
    for (int i = 0; i < nCrystals; ++i) {
        CalGeEnergy[i] = new TH1F(Form("CalGeEnergy%02d", i),Form("Cal_e - %02d", i), 4000, 0, 4000);
        BGOGeTimeDiff[i] = new TH1F(Form("BGOGeTimeDiff%02d", i),Form("BGOGeTimeDiff_Plot%02d", i),600,-3000,3000);
        CalGeEnergy[i]->SetDirectory(nullptr);//comment this out if you want to use this histogram
    }
    // Create addback gamma energy histograms - ONE PER CLOVER
    TH1F* AddBackGeEnergy[nClovers];
    for (int i = 0; i < nClovers; ++i) {
        AddBackGeEnergy[i] = new TH1F(Form("AddBackGeEnergy%02d", i),Form("AddBack_e - %02d", i), 4000, 0, 4000);
    }
    //Create summed clover energy spectra (just summed not addback)
    TH1F* CompSupCloverSumEnergy[nClovers];
    for (int i = 0; i < nClovers; ++i) {
        CompSupCloverSumEnergy[i] = new TH1F(Form("SummedCloverCalGeEnergy%02d", i),Form("Summed_e - %02d", i), 4000, 0, 4000);
        CompSupCloverSumEnergy[i]->SetDirectory(nullptr);//comment this out if you want to use this histogram
    }
    //Create Ge Time difference plot in each clover
    TH1F* GeTimeDiff[nCrystalPairs];
    for (int i = 0; i<nCrystalPairs; i++) {
        GeTimeDiff[i] = new TH1F(Form("GeTimeDiff%02d", i),Form("GeTimeDiff_Plot%02d", i),600,-3000,3000);
    }

    //Create Labr Time difference plot
    TH1F* LabrTimeDiff[nLabrPairs];
    for (int i = 0; i<nLabrPairs; i++) {
        LabrTimeDiff[i] = new TH1F(Form("LabrTimeDiff%02d", i),Form("LabrTimeDiff_Plot%02d", i),3000,0,3000);
        LabrTimeDiff[i]->SetOption("HIST");
    }
    
    //Create Labr Time difference plot vs labr energy
    TH2F* LabrTimeDiffvLabrEnergy[nLabrPairs];
    for (int i = 0; i<nLabrPairs; i++) {
        LabrTimeDiffvLabrEnergy[i] = new TH2F(Form("LabrTimeDiffvLabrEnergy%02d", i),Form("LabrTimeDiffvLabrEnergy_Plot%02d", i),200,-100,100,3000,0,3000);
        //LabrTimeDiffvLabrEnergy[i]->SetOption("COLZ");
    }
    
    //Create a total time diff plot Ge
    TH1F* TotalGeTimeDiff = new TH1F("TotalGeTimeDiff","TotalGeTimeDiff_plot",600,-3000,3000);

    //Create a total time diff plot Ge
    TH1F* TotalLaBr3TimeDiff = new TH1F("TotalLaBr3TimeDiff","TotalLaBr3TimeDiff_plot",48000,-1000,1000);
    TotalLaBr3TimeDiff->SetDirectory(nullptr);//comment this out if you want to use this histogram

    //Create LaBr3 energy hists
    TH1F* LaBr3Energy[nLaBr3];
    for (int i = 0; i < nLaBr3; ++i) {
        LaBr3Energy[i] = new TH1F(Form("LaBr3Energy%02d", i),Form("e - %02d", i), 6000, 0, 10000);
        //LaBr3Energy[i]->SetDirectory(nullptr);//comment this out if you want to use this histogram
        //LaBr3Energy[i]->SetOption("HIST");
    }
    
    //Create Raw BEGe energy hists
    TH1F* rawBEGe[nBEGe];
    for (int i = 0; i < nBEGe; ++i) {
        rawBEGe[i] = new TH1F(Form("rawBEGe%02d", i),Form("e - %02d", i), 4000, 0, 16000);
        rawBEGe[i]->SetDirectory(nullptr);//comment this out if you want to use this histogram
    }
    //Create Raw BEGe energy hists
    TH1F* calBEGe[nBEGe];
    for (int i = 0; i < nBEGe; ++i) {
        calBEGe[i] = new TH1F(Form("calBEGe%02d", i),Form("e - %02d", i), 4000, 0, 4000);
        calBEGe[i]->SetDirectory(nullptr);//comment this out if you want to use this histogram
    }
    //Create the PID hists
    //std::cout<<"Before PIDs"<<std::endl;
    //TH2F* hPID[NGAGG];
    //TH2F* hPID_alpha[NGAGG];
    //for (int i = 0; i<NGAGG; i++){
       // hPID[i] = new TH2F(Form("hPID%02d",i), Form("PID - %02d",i), PIDnBins, 0, PIDxSize, PIDnBins, 0, PIDySize);
       // hPID_alpha[i] = new TH2F(Form("hPID_alpha%02d",i), Form("PID_alpha - %02d",i), PIDnBins, 0, PIDxSize, PIDnBins, 0, PIDySize);
    //}
    //Create Alpha Cut array
    TCutG* alphaCuts[NGAGG] = {nullptr};

    //Create a hist for gamma gamma
    TH2F* CompSupGammaGamma = new TH2F("CompSupGammaGamma","CompSupGammaGamma_plot",1750,0,4000,1750,0,4000);
    CompSupGammaGamma->SetDirectory(nullptr);//comment this out if you want to use this histogram

    //Create hist for time gated gamma gamma
    TH2F* TotalGammaGamma = new TH2F("TotalGammaGamma","TotalGammaGamma_plot",1750,0,4000,1750,0,4000);
    TotalGammaGamma->SetDirectory(nullptr);//comment this out if you want to use this histogram

    //Create hist for addback gamma gamma
    TH2F* AddBackGammaGamma = new TH2F("AddBackGammaGamma","AddBackGammaGamma_plot",1750,0,4000,1750,0,4000);

    //Create hist for alphaCut addback gamma gamm
    TH2F* alphaCutAddBackGammaGamma = new TH2F("alphaCutAddBackGammaGamma","alphaCutAddBackGammaGamma_plot",1750,0,4000,1750,0,4000);

    //Create a hist to hold the 500MHz cfd source bits
    TH1F* forcedTrigger500MHzBoard = new TH1F("forcedTrigger500MHzBoard","forcedTrigger500MHzBoard",11,-1,10);

    //Create hists for traces from each detector
    TH2F* Traces[nLaBr3];
    for (int i = 0; i < nLaBr3; ++i) {
        Traces[i] = new TH2F(Form("traces%02d", i),Form("traces - %02d", i), 600, 0, 300, 10000, 0, 10000);
    }




    //This reads in the alphaCuts root file
    readAlphaCutsFile(alphaCutsFile, alphaCuts);


    // Fill gamma energy histograms
    Long64_t nentries = tree->GetEntries();
    for (Long64_t i = 0; i < nentries; ++i) {
        tree->GetEntry(i);
        //This fills the Raw (and wriggling!) germanium energy histograms
        makeWrigglingGerm(tree,kMaxMulti,multi,detNum,e,RawGeEnergy);

        //This fills the calibrated Ge energy
        makeCalGeEnergy(tree,kMaxMulti,multi,detNum,e,CalGeEnergy);

        //This fills the AddBack calibrated Ge Energy
        makeAddBackGeEnergy(tree,kMaxMulti,multi,detNum,e,cloverNum,AddBackGeEnergy, includeDiagonals);

        //This fills the PIDs
        //makePID(tree,kMaxMulti,NGAGG,multi,detNum,qdc,hPID,trinCoin);

        //This fills the BGOs
        //makeBGO(tree,kMaxMulti,multi,detNum,e,BGOEnergy);

        //This fills the LaBr3 energy
        makeLaBr3(tree, kMaxMulti, multi, detNum,e,LaBr3Energy,qdc);

        //This fills the Raw BEGe
        makeRawBEGe(tree, kMaxMulti, multi, detNum,e,rawBEGe);

        //This fills the calibrated BEGe
        //makeCalBEGe(tree, kMaxMulti, multi, detNum,e,calBEGe);

        //this fills the Compton suppressed Ge Energy hists
        makeCompSupCalGeEnergy(tree,kMaxMulti,multi,detNum,e_t,e,CompSupCalGeEnergy,BGOGeTimeGate_ns);

        //This fills a compton suppressed gamma gamma matrix across all clovers]
        //makeCompSupCalGammaGamma(tree,kMaxMulti,multi,detNum,e_t,e,CompSupGammaGamma,cloverNum,BGOGeTimeGate_ns);

        //This fills a compton suppressed gamma gamma matrix across all clovers with a Ge time cut cats
        //makeCalGammaGamma(tree,kMaxMulti,multi,detNum,e_t,e,TotalGammaGamma,BGOGeTimeGate_ns);

        //This fills an addback'ed gamma gamma
        //makeAddBackGammaGamma(tree,kMaxMulti,multi,detNum,e_t,e,AddBackGammaGamma,BGOGeTimeGate_ns,cloverNum,pileup);

        //This fills a time difference plot per clover
        makeGeTimeDiff(tree,kMaxMulti,multi,detNum,e_t,e,GeTimeDiff,cloverNum,pileup);

        //This filles a time difference plot per labr pair
        makeLabrTimeDiff(tree, kMaxMulti, multi, detNum, LabrTimeDiff, pileup, cfd, qdc, e_t);

        //This filles a time difference plot per labr pair vs the first labr energy of the pair found in event.
        makeLabrTimeDiffvLabrEnergy(tree, kMaxMulti, multi, detNum, LabrTimeDiffvLabrEnergy, pileup, cfd, qdc, e_t);

        //This fills a Ge time difference plot regardless of detector/clover. This should just be the sum of the above.
        makeTotalGeTimeDiff(tree,kMaxMulti,multi,detNum,e_t,e,TotalGeTimeDiff,cloverNum);

        //This fills a LaBr3 time difference plot regardless of detector.
        makeTotalLaBr3TimeDiff(tree,kMaxMulti,multi,detNum,e_t,e,TotalLaBr3TimeDiff,cloverNum,cfd);

        //This files a BGO-Ge time diff plot regardless of detector
        makeBGOGeTimeDiff(tree,kMaxMulti,multi,detNum,e_t,e,BGOGeTimeDiff,cloverNum);

        //This fills the alphacut addback gamma gamma mat
        //makeAlphaCutAddBackGammaGamma(tree,kMaxMulti,multi,detNum,qdc,e_t,e,alphaCutAddBackGammaGamma,BGOGeTimeGate_ns,cloverNum,pileup,alphaCuts);

        //This fills the 500MHz board cfd source bits into a hist. Diagnostic to check if code is throwing some out.
        makeCFDSourceBits(multi, detNum, forcedTrigger500MHzBoard, cfd);

        //This fills the traces histograms
        makeTraces(multi, detNum, trace, Traces);



        //Progress Bar
        if (i % (nentries / 100) == 0) {
            float progress = (float)i / (float)nentries;
            int pos = barWidth * progress;
        
            std::cout << "\r[";
            for (int k = 0; k < barWidth; ++k) {
                if (k < pos)
                    std::cout << "\033[42m \033[0m";         // Green fill
                else
                    std::cout << "\033[41m \033[0m";         // Red space
            }
            std::cout << "] " << int(progress * 100.0) << "% "
                      << "\033[31m" << spinner[spinIdx++ % 4] << "\033[0m" << std::flush;
        }  
    }

    //This makes a singles sum energy comp sup per clover
    //makeCompSupCloverSumEnergy(CompSupCalGeEnergy, CompSupCloverSumEnergy);

    //******************************************************************************Cuts Being Applied here!!!!***************************************************** */

    //This fills the alpha cut PIDs
    //makeAlphaCutPID(hPID,alphaCuts, hPID_alpha);

    

    //Progress bar
    std::cout << "\r[";
    for (int k = 0; k < barWidth; ++k)
        std::cout << "\033[42m \033[0m";  // Full green bar
    std::cout << "] 100% \033[32m✓\033[0m\n" << std::endl;  // Red checkmark
    std::cout <<"\033[31mPlease wait for histograms to fill!!!!\033[0m\n"<<std::endl;
    
    
    // Save histograms to output file
    std::string outFileName = std::string(outputDirectory) + std::string(outputFilename);
    TFile *outputFile = new TFile(outFileName.c_str(), "RECREATE");
    for (int i = 0; i < nCrystals; ++i) {
        RawGeEnergy[i]->Write();
        CompSupCalGeEnergy[i]->Write();
        //CalGeEnergy[i]->Write();
    }
    //for (int i = 0; i < nClovers; ++i) {
        //CompSupCloverSumEnergy[i]->Write();
    //}
    //for (int i = 0; i < nClovers; ++i) {
        //BGOEnergy[i]->Write();
   //}
    for (int i = 0; i < nClovers; ++i) {
        AddBackGeEnergy[i]->Write();
    }
    //for (int i = 0; i < nZeroDeg; ++i) {
        //ZeroDegEnergy[i]->Write();
    //}
    //for (int i = 0; i < NGAGG; ++i) {
        //hPID[i]->Write();
    //}
    //for (int i = 0; i < NGAGG; ++i) {
        //hPID_alpha[i]->Write();
    //}
    for (int i = 0; i < nLaBr3; ++i) {
        LaBr3Energy[i]->Write();
        //std::cout<<"I WROTE THE LABR"<<std::endl;
    }
    for (int i = 0; i < nBEGe; ++i) {
        rawBEGe[i]->Write();
        //calBEGe[i]->Write();
    }
    for(int i = 0; i<nCrystalPairs;i++){
        GeTimeDiff[i]->Write();
    }
    for(int i = 0; i<nLabrPairs;i++){
        LabrTimeDiff[i]->Write();
        LabrTimeDiffvLabrEnergy[i]->Write();
    }
    for(int i = 0; i<nCrystals; i++){
        BGOGeTimeDiff[i]->Write();
    }
    for(int i =0; i<nLaBr3; i++){
        Traces[i]->Write();
    }
    
    TotalGeTimeDiff->Write();
    TotalLaBr3TimeDiff->Write();
    //CompSupGammaGamma->Write();
    //TotalGammaGamma->Write();
    AddBackGammaGamma->Write();
    //alphaCutAddBackGammaGamma->Write();
    forcedTrigger500MHzBoard->Write();
    outputFile->Close();

    // Cleanup
    delete outputFile;

    std::cout << "Histograms Saved!\n";
}

int main(int argc, char** argv){
    if (argc<12){
        std::cerr<<"Usage: "<<argv[0]
                            <<" <input Data File>"
                            <<" <calFile.dat>"
                            <<" <gePromptGatesFile.dat>"
                            <<" <output Directory>"
                            <<" <output Filename (include extension)>"
                            <<" <enforce trinity A and B coin? (1 = yes 0 = No)>"
                            << " <include diagonals in addback? (1 = yes 0 = no)>"
                            <<" <alphaCuts File name (inlcude .root)>"
                            <<" <geTimeingPairFile (include extension)>"
                            <<" <BGO-Ge time gate>"
                            <<" <labrTimingPairFile (include extension)>"
                            <<std::endl;
        return 1;
    }
    bool trinCoin = 0;
    bool includeDiagonals = 0;
    const char* dataFile = argv[1];
    const char* calFile = argv[2];
    const char* gePromptGatesFile = argv[3];
    const char* outDir = argv[4];
    const char* outFile = argv[5];
    const char* trinCoinChar = argv[6];
    const char* includeDiagonalsChar = argv[7];
    const char* alphaCutsFile = argv[8];
    const char* geTimingPairFile = argv[9];
    const char* BGOGeTimeCut_ns = argv[10];
    const char* labrTimingPairFile = argv[11];

    if (std::string(trinCoinChar) == "1"){
        trinCoin = true;
        std::cout<<"Trinity A and B coincidence \033[32mENFORCED\033[0m"<<std::endl;
    }else if (std::string(trinCoinChar) == "0"){
        trinCoin = false;
        std::cout<<"Trinity A and B coincidence \033[31mNOT\033[0m ENFORCED"<<std::endl;
    }
    if (std::string(includeDiagonalsChar) == "1"){
        includeDiagonals = true;
        std::cout<<"Diagonls are \033[32mINCLUDED\033[0m in addback"<<std::endl;
    }else if (std::string(includeDiagonalsChar) == "0"){
        includeDiagonals = false;
        std::cout<<"Diagonals are \033[31mNOT\033[0m INCLUDED in addback"<<std::endl;
    }
    std::cout<<"Applying a \033[1;36m"<<BGOGeTimeCut_ns<<"\033[0m nanosecond BGO-Ge time gate"<<std::endl;
    Int_t BGOGeTimeCut_ns_int = std::stoi(BGOGeTimeCut_ns);
    
    //std::cout<<"I am right before the Process_Events call in main"<<std::endl;
    Process_Events(dataFile, calFile, gePromptGatesFile, outDir, outFile, trinCoin, includeDiagonals, alphaCutsFile, geTimingPairFile, BGOGeTimeCut_ns_int, labrTimingPairFile);
    return 0;
}