rtang/XIAEventBuilder
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
XIAEventBuilder/JakeStuff/Jericho.cpp

1568 lines
69 KiB
C++
Raw Permalink Blame History

/******************************************
* 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;
}
Reference in New Issue View Git Blame Copy Permalink