//Macro for subtracting gated spectra for coincidence analysis. It first defines new 1D hist and then subtracts the background spectra from the gated spectrum. This will then subtract out the gated peak to then analyze the remaining peaks in coincidence within the observed matrix.

//x means that the line works as it shoulds - work in progress (4/12/23)

#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TCutG.h>
#include <TFile.h>
#include "TH2.h"
#include "TH1.h"
#include "TRandom.h"
#include "iostream"
#include <TMath.h>
#include "TObjArray.h"
#include "TF1.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"
#include "Riostream.h"

TObjArray List1;
TObjArray List2;
TObjArray List3;

void bgs(){

//Place in multiple histograms into the same canvas (usage of pads):

//10(width),10(width), 800(height), 800(height)
//TCanvas *c1 = new TCanvas("c1", "Gated Spectra", 10,10, 800, 800);
//i.e. pad = new TPad(name,name,x1,y1,x2,y2,0);
//TPad *pad1 = new TPad("pad1","top pad",0.2,0.6,0.8,0.99);
//pad1->Draw();
//pad1->cd();
//TPad *pad2 = new TPad("pad2","bottom pad",0,0,1,0.5);
//pad2->Draw();
//pad2->cd();
//TCanvas *c2 = new TCanvas("c2","Avg_Bg Spectra",10,10, 800, 800);
//TCanvas *c3 = new TCanvas("c3", "Gated Minus Avg_Bg", 10,10, 800, 800);

//func = G - 0.5(background_gate + background_gate_2) -> function used to subtract out the gated peak of interest (POI).

//Retrieve the histogram:

TFile *file = new TFile("out_gate_1226_1526_1026_1126_1626_1726.root");

//Placing multiple histograms in the same canvas:

TCanvas *c0 = new TCanvas("c0", "64Ni Coincidence Analysis", 10,10, 800, 800);
//(columns,rows) -> c0
c0->Divide(2,2);
c0->cd(1);

//TH1D *his = file->Get("histoFit");
//4096
//16000

TH1D * gated=new TH1D("gate","gate",16000,0,16000);
TH1D * avg_bg=new TH1D("avg_bg","avg_bg",4096,0,4096);
TH1D * gated_bgs=new TH1D("gated_bgs","gated_bgs",16000,0,16000);
TH1D * gated_tbgs=new TH1D("gated_tbgs","gated_tbgs",16000,0,16000);

//Adding the gated spectrum first to the empty canvas (pad1):

gated->Add(gate,1);
gated->Draw();

//Placing the background regions into the canvas then subtract the gated_spectrum from the avg_bg (scale factor determined in hdtv):

//Adding the background spectra (pad2):

c0->cd(2);
avg_bg->Add(background_gate,0.5);
avg_bg->Add(background_gate_2,0.5);
avg_bg->Draw();

//Subtracting gate from avg_bg into pad3:

c0->cd(3);
gated_bgs->Add(gate,1);
Double_t scale_1 = 1.651;//This is followed by scaling the avg_bg. 
avg_bg->Scale(scale_1);
gated_bgs->Add(avg_bg,-1);
gated_bgs->Draw();

//Subtracting the gated_spectrum from the avg_bg again to remove the gated peak. Applying a second scale to align the gated spectrum with the average background spectrum (using the pair production - 511 MeV peak [scale factor determined in hdtv] - pad4):
c0->cd(4);
gated_tbgs->Add(gated_bgs,1);
Double_t scale_2 = 0.283;
avg_bg->Scale(scale_2);
gated_tbgs->Add(avg_bg,1);
gated_tbgs->Draw();

}

//3 parameters for Peak of Interest (POI) and Background. Setup the TF1 functions:

//TF1 *f1 = new TF1("f1","p0*exp(-0.5*((x-p1)/p2)^2)",0,16000);NOTE the range! (edit parameters) 
//TF1 *func = new TF1("func","2.84626e+03*exp(-0.5*((x-1.32539e+03)/4.16607e+01)^2)",0,16000);

//Fitting the gaussian parameters for the gate and background:

//TF1 *fGate = new TF1("gate","gaus",1250.,30.);
//TF1 *fAvgBackground = new TF1("avg_bg","pol0+expo(1)",1250.,.);
//TF1 *fBackground2 = new TF1("background_gate2","pol0+expo(1)",0.,30.); 
//TF1 *fSpectrum = new TF1("tbs","gaus+pol0(3)+expo(4)",0.,30.);

//6 parameters for POI + Background. "Strength","Mean","Sigma"-> provided from peak fitting in ROOT.
//fSpectrum->SetParNames("Strength","Mean","Sigma","Back1","Back2","Back3");
//fSpectrum->SetParameters(2.84626e+03, 1.32539e+03, 4.16607e+01, 16000, 0, 16000);

//Fit the spectrum:

//his->Fit("fSpectrum”, ””, ””, 0., 30.);
//Double_t param[6];

//Set all TF1 functions to fitted parameters:

//fSpectrum->GetParameters(param);
//fSignal->SetParameters(&param[0]);
//fBackground->SetParameters(&param[3]);

//SUbtracts the background:

//TH1D *hisSignal = new TH1D(*his);
//hisSignal->Sumw2();
//hisSignal->Add(fBackground, -1);

//Plot the spectra and functions:

//his->Draw("e"); hisSignal->Draw(“SAME”);
//fSignal->Draw("same"); fBackground->Draw("same");