rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

modified: GainMatchSX3.C

Browse Source 
modified:   GainMatchSX3Front.C
    changes made to the GainMatchSX3.C and GainMatchSX3Front.C files to include reduced Chisquared and fixed uncertainties for the gain matching process.
This commit is contained in:
Vignesh Sitaraman 2025-07-24 16:12:46 -04:00
parent 4fc05ea338
commit 3d0d176f5a
2 changed files with 225 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

189
GainMatchSX3.C
View File

@ -12,7 +12,7 @@
#include <algorithm>
#include <TProfile.h>
#include "Armory/ClassSX3.h"
#include <TGraphErrors.h>
#include "TVector3.h"
TH2F *hSX3FvsB;
@ -26,7 +26,17 @@ int padID = 0;
SX3 sx3_contr;
TCutG *cut;
TCutG *cut1;
std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
std::map<std::tuple<int, int, int, int>, int> comboCounts;
const int MAX_DET = 24;
const int MAX_UP = 4;
const int MAX_DOWN = 4;
const int MAX_BK = 4;
double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
void GainMatchSX3::Begin(TTree * /*tree*/)
{
@ -56,6 +66,33 @@ void GainMatchSX3::Begin(TTree * /*tree*/)
return;
}
cut->SetName("sx3cut"); // Ensure the cut has the correct name
// Load the TCutG object
TFile *cutFile1 = TFile::Open("UvD.root");
bool cut1Loaded = (cut1 != nullptr);
cut1 = dynamic_cast<TCutG *>(cutFile1->Get("UvD"));
if (!cut1)
{
std::cerr << "Error: Could not find TCutG named 'UvD' in UvD.root" << std::endl;
return;
}
cut1->SetName("UvD");
std::string filename = "sx3_GainMatchfront.txt";
std::ifstream infile(filename);
if (!infile.is_open())
{
std::cerr << "Error opening " << filename << "!" << std::endl;
return;
}
int id, bk, u, d;
double gain;
while (infile >> id >> bk >> u >> d >> gain)
{
frontGain[id][bk][u][d] = gain;
frontGainValid[id][bk][u][d] = true;
}
}
Bool_t GainMatchSX3::Process(Long64_t entry)
@ -135,29 +172,38 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
for (size_t i = 0; i < sx3ID.size(); i++)
{
int index = sx3ID[i].second;
// Check the channel number and assign it to the appropriate channel type
if (sx3.ch[index] < 8)
if (sx3.e[i] > 100)
{
if (sx3.ch[index] % 2 == 0)
int index = sx3ID[i].second;
// Check the channel number and assign it to the appropriate channel type
if (sx3.ch[index] < 8)
{
sx3ChDn = sx3.ch[index];
sx3EDn = sx3.e[index];
if (sx3.ch[index] % 2 == 0)
{
sx3ChDn = sx3.ch[index];
sx3EDn = sx3.e[index];
}
else
{
sx3ChUp = sx3.ch[index];
sx3EUp = sx3.e[index];
}
}
else
{
sx3ChUp = sx3.ch[index];
sx3EUp = sx3.e[index];
sx3ChBk = sx3.ch[index] - 8;
// if (sx3ChBk == 2)
// printf("Found back channel Det %d Back %d \n", sx3.id[index], sx3ChBk);
sx3EBk = sx3.e[index];
}
}
else
{
sx3ChBk = sx3.ch[index] - 8;
// if (sx3ChBk == 2)
// printf("Found back channel Det %d Back %d \n", sx3.id[index], sx3ChBk);
sx3EBk = sx3.e[index];
}
}
for (int i = 0; i < sx3.multi; i++)
{
auto key = std::make_tuple(sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn);
comboCounts[key]++;
}
// If we have a valid front and back channel, fill the histograms
hSX3->Fill(sx3ChDn, sx3ChBk);
hSX3->Fill(sx3ChUp, sx3ChBk);
@ -167,8 +213,13 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
for (int i = 0; i < sx3.multi; i++)
{
if (sx3.id[i] == 3)
// if (sx3.id[i] == 4)
{
auto key = std::make_tuple(sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn);
// Only continue if this combo has enough entries
if (comboCounts[key] < 100 || sx3EBk < 100 || sx3EUp < 100 || sx3EDn < 100)
continue;
// Fill the histogram for the front vs back with gain correction
hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
// Fill the index vs energy histogram
@ -189,10 +240,13 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
// if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk))
// if (sx3.id[i] < 24 && sx3ChUp < 4 && sx3ChBk < 4 && std::isfinite(sx3EUp) && std::isfinite(sx3EDn) && std::isfinite(sx3EBk))
if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk))// && cut1 && cut1->IsInside(sx3EUp / sx3EBk, sx3EDn / sx3EBk))
{
// Accumulate data for gain matching
// if (frontGainValid[sx3.id[i]][sx3ChBk][sx3ChUp][sx3ChDn])
// {
// sx3EUp *= frontGain[sx3.id[i]][sx3ChBk][sx3ChUp][sx3ChDn];
// }
dataPoints[{sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(sx3EBk, sx3EUp, sx3EDn);
}
}
@ -245,47 +299,88 @@ void GainMatchSX3::Terminate()
{
double eUp, eDn, eBk;
std::tie(eBk, eUp, eDn) = pr;
if ((eBk < 100) || (eUp < 100) || (eDn < 100))
continue; // Skip if any energy is less than 100
bkE.push_back(eBk);
udE.push_back(eUp + eDn);
}
// Fill the TGraph with bkE and udE
TGraph g(bkE.size(), bkE.data(), udE.data());
// TGraph g(bkE.size(), bkE.data(), udE.data());
// Fit the graph to a linear function
TF1 f("f", "[0]*x", 0, 16000);
g.Fit(&f, "QNR");
if (TMath::Abs(f.GetParameter(0) - 1) > 1)
if (bkE.size() < 5)
continue; // Ensure we have enough points for fitting
// TF1 f("f", "[0]*x", 0, 16000);
// g.Fit(&f, "NR");
// if (TMath::Abs(f.GetParameter(0) - 1) > 3.0)
// continue;
const double fixedError = 20.0; // in ADC channels
std::vector<double> xVals, yVals, exVals, eyVals;
// Build data with fixed error
for (size_t i = 0; i < udE.size(); ++i)
{
continue; // Skip this fit if the slope is too far from 1
double x = udE[i]; // front energy
double y = bkE[i]; // back energy
xVals.push_back(x);
yVals.push_back(y);
// exVals.push_back(fixedError); // error in front energy
eyVals.push_back(fixedError); // error in back energy
}
// Build TGraphErrors with errors
TGraphErrors g(xVals.size(), xVals.data(), yVals.data(), exVals.data(), eyVals.data());
TF1 f("f", "[0]*x", 0, 16000);
// f.SetParameter(0, 1.0); // Initial guess
// Interactive canvas
TCanvas *c = new TCanvas(Form("c_%d_%d_%d_%d", id, bk, u, d), "Fit", 800, 600);
g.SetTitle(Form("Detector %d: U%d D%d B%d", id, u, d, bk));
g.SetMarkerStyle(20);
g.SetMarkerColor(kBlue);
g.Draw("AP");
g.Fit(&f, "Q"); // Quiet fit
double chi2 = f.GetChisquare();
int ndf = f.GetNDF();
double reducedChi2 = (ndf != 0) ? chi2 / ndf : -1;
std::cout << Form("Det%d U%d D%d B%d → Gain: %.4f | χ²/ndf = %.2f/%d = %.2f",
id, u, d, bk, f.GetParameter(0), chi2, ndf, reducedChi2)
<< std::endl;
// Show canvas and wait for user to continue
c->Update();
gPad->WaitPrimitive();
if (TMath::Abs(f.GetParameter(0) - 1) > 3.0)
{
delete c; // Clean up unused canvas
continue;
}
gainArray[id][bk][u][d] = f.GetParameter(0);
gainValid[id][bk][u][d] = true;
}
// Output results
for (int id = 0; id < MAX_DET; ++id)
{
for (int bk = 0; bk < MAX_BK; ++bk)
// Check if the gain is valid for this detector, back, up, and down
if (gainValid[id][bk][u][d])
{
for (int u = 0; u < MAX_UP; ++u)
if (TMath::Abs(gainArray[id][u][d][bk] - 1) < 0.3)
{
for (int d = 0; d < MAX_DOWN; ++d)
{
// Check if the gain is valid for this detector, back, up, and down
if (gainValid[id][bk][u][d])
{
outFile << id << " " << bk << " " << u << " " << d << " " << gainArray[id][u][d][bk] << std::endl;
if (TMath::Abs(gainArray[id][u][d][bk] - 1) < 0.3)
{
printf("Gain match Det%d Up%dDn%d Back%d → %.4f \n", id, u, d, bk, gainArray[id][u][d][bk]);
}
else if (gainArray[id][u][d][bk]!=0)
{
std::cerr << "Warning: Gain value out of range for Det " << id << " Up " << u << " Dn " << d << " Back " << bk << ": "
<< gainArray[id][u][d][bk] << std::endl;
}
}
}
printf("Gain match Det%d Up%dDn%d Back%d → %.4f \n", id, u, d, bk, gainArray[id][u][d][bk]);
outFile << id << " " << bk << " " << u << " " << d << " " << gainArray[id][u][d][bk] << std::endl;
}
else if (gainArray[id][u][d][bk] != 0)
{
std::cerr << "Warning: Gain value out of range for Det " << id << " Up " << u << " Dn " << d << " Back " << bk << ": "
<< gainArray[id][u][d][bk] << std::endl;
}
}
}

117
GainMatchSX3Front.C
View File

@ -91,7 +91,10 @@ void GainMatchSX3Front::Begin(TTree * /*tree*/)
while (infile >> id >> bk >> u >> d >> gain)
{
backGain[id][bk][u][d] = gain;
backGainValid[id][bk][u][d] = true;
if(backGain[id][bk][u][d] > 0)
backGainValid[id][bk][u][d] = true;
else
backGainValid[id][bk][u][d] = false;
}
infile.close();
@ -116,7 +119,7 @@ Bool_t GainMatchSX3Front::Process(Long64_t entry)
for (int j = i + 1; j < sx3.multi; j++)
{
if (sx3.id[i] == 3)
// if (sx3.id[i] == 3)
hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
}
if (sx3.e[i] > 100)
@ -196,7 +199,7 @@ Bool_t GainMatchSX3Front::Process(Long64_t entry)
for (int i = 0; i < sx3.multi; i++)
{
if (sx3.id[i] == 3 && sx3.e[i] > 100)
if ( sx3.e[i] > 100)// && sx3.id[i] == 4)
{
// back gain correction
@ -220,8 +223,8 @@ Bool_t GainMatchSX3Front::Process(Long64_t entry)
hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk) &&
cut1 && cut1->IsInside(sx3EUp / sx3EBk, sx3EDn / sx3EBk))
if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk)
&& cut1 && cut1->IsInside(sx3EUp / sx3EBk, sx3EDn / sx3EBk))
{
if (backGainValid[sx3.id[i]][sx3ChBk][sx3ChUp][sx3ChDn])
@ -231,11 +234,6 @@ Bool_t GainMatchSX3Front::Process(Long64_t entry)
// Accumulate data for gain matching
dataPoints[{sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(sx3EBk, sx3EUp, sx3EDn);
}
// if (sx3.id[i] < 24 && sx3ChUp < 4 && sx3ChBk < 4 && std::isfinite(sx3EUp) && std::isfinite(sx3EDn) && std::isfinite(sx3EBk))
{
// Accumulate data for gain matching
dataPoints[{sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(sx3EBk, sx3EUp, sx3EDn);
}
}
}
}
@ -268,20 +266,71 @@ void GainMatchSX3Front::Terminate()
if (pts.size() < 5)
continue;
std::vector<double> udE, corrBkE;
std::vector<double> uE, dE, udE, corrBkE;
for (const auto &pr : pts)
{
double eBkCorr, eUp, eDn;
std::tie(eBkCorr, eUp, eDn) = pr;
if( (eBkCorr < 100) || (eUp <100) || (eDn < 100))
continue; // Skip if any energy is zero
uE.push_back(eUp / eBkCorr);
dE.push_back(eDn / eBkCorr);
udE.push_back(eUp + eDn);
corrBkE.push_back(eBkCorr);
hUvD->Fill(eUp / eBkCorr, eDn / eBkCorr);
}
if( uE.size() < 5 || dE.size() < 5 || corrBkE.size() < 5)
continue; // Ensure we have enough points for fitting
// TGraph g(udE.size(), udE.data(), corrBkE.data());
// TF1 f("f", "[0]*x", 0, 20000);
// f.SetParameter(0, 1.0); // Initial guess for the gain
// g.Fit(&f, "R");
const double fixedError = 20.0; // in ADC channels
std::vector<double> xVals, yVals, exVals, eyVals;
// Build data with fixed error
for (size_t i = 0; i < udE.size(); ++i)
{
double x = udE[i]; // front energy
double y = corrBkE[i]; // back energy
xVals.push_back(x);
yVals.push_back(y);
// exVals.push_back(fixedError); // error in front energy
eyVals.push_back(fixedError); // error in back energy
}
// Build TGraphErrors with errors
TGraphErrors g(xVals.size(), xVals.data(), yVals.data(), exVals.data(), eyVals.data());
TF1 f("f", "[0]*x", 0, 16000);
// f.SetParameter(0, 1.0); // Initial guess
// Interactive canvas
TCanvas *c = new TCanvas(Form("c_%d_%d_%d_%d", id, bk, u, d), "Fit", 800, 600);
g.SetTitle(Form("Detector %d: U%d D%d B%d", id, u, d, bk));
g.SetMarkerStyle(20);
g.SetMarkerColor(kBlue);
g.Draw("AP");
g.Fit(&f, "Q"); // Quiet fit
double chi2 = f.GetChisquare();
int ndf = f.GetNDF();
double reducedChi2 = (ndf != 0) ? chi2 / ndf : -1;
std::cout << Form("Det%d U%d D%d B%d → Gain: %.4f | χ²/ndf = %.2f/%d = %.2f",
id, u, d, bk, f.GetParameter(0), chi2, ndf, reducedChi2)
<< std::endl;
// Show canvas and wait for user to continue
c->Update();
gPad->WaitPrimitive();
TGraph g(udE.size(), udE.data(), corrBkE.data());
TF1 f("f", "[0]*x", 0, 40000);
g.Fit(&f, "QNR");
frontGain[id][bk][u][d] = f.GetParameter(0);
frontGainValid[id][bk][u][d] = true;
@ -294,7 +343,7 @@ void GainMatchSX3Front::Terminate()
std::cout << "Gain matching complete." << std::endl;
// === Stage 3: Create corrected histogram ===
TH2F *hCorrectedFvB = new TH2F("hCorrectedFvB", "Corrected;Corrected Front Sum;Corrected Back", 800, 0, 16000, 800, 0, 16000);
TH2F *hCorrectedFvB = new TH2F("hCorrectedFvB", "Corrected;Corrected Front Sum;Corrected Back", 800, 0, 8000, 800, 0, 8000);
TH2F *hCorrectedUvD = new TH2F("hCorrectedUvD", "Corrected UvD; UvD Up; UvD Down", 600, 0, 1, 600, 0, 1);
for (const auto &kv : dataPoints)
@ -308,31 +357,31 @@ void GainMatchSX3Front::Terminate()
double eBk, eUp, eDn;
std::tie(eBk, eUp, eDn) = pr;
double corrUp = eUp * front;
double corrDn = eDn * front;
// double corrDn = eDn * front;
hCorrectedFvB->Fill(corrUp + corrDn, eBk);
hCorrectedUvD->Fill(corrUp / eBk, corrDn / eBk);
hCorrectedFvB->Fill(corrUp + eDn, eBk);
hCorrectedUvD->Fill(corrUp / eBk, eDn / eBk);
}
}
// === Final canvas ===
gStyle->SetOptStat(1110);
TCanvas *c1 = new TCanvas("c1", "Gain Correction Results", 1200, 600);
c1->Divide(2, 1);
// // === Final canvas ===
// gStyle->SetOptStat(1110);
// TCanvas *c1 = new TCanvas("c1", "Gain Correction Results", 1200, 600);
// c1->Divide(2, 1);
c1->cd(1);
hSX3FvsB_g->SetTitle("Before Correction (Gated)");
hSX3FvsB_g->GetXaxis()->SetTitle("Measured Front Sum (E_Up + E_Dn)");
hSX3FvsB_g->GetYaxis()->SetTitle("Measured Back E");
hSX3FvsB_g->Draw("colz");
// c1->cd(1);
// hSX3FvsB_g->SetTitle("Before Correction (Gated)");
// hSX3FvsB_g->GetXaxis()->SetTitle("Measured Front Sum (E_Up + E_Dn)");
// hSX3FvsB_g->GetYaxis()->SetTitle("Measured Back E");
// hSX3FvsB_g->Draw("colz");
c1->cd(2);
hCorrectedFvB->SetTitle("After Correction");
hCorrectedFvB->Draw("colz");
TF1 *diag = new TF1("diag", "x", 0, 40000);
diag->SetLineColor(kRed);
diag->SetLineWidth(2);
diag->Draw("same");
// c1->cd(2);
// hCorrectedFvB->SetTitle("After Correction");
// hCorrectedFvB->Draw("colz");
// TF1 *diag = new TF1("diag", "x", 0, 40000);
// diag->SetLineColor(kRed);
// diag->SetLineWidth(2);
// diag->Draw("same");
std::cout << "Terminate() completed successfully." << std::endl;
}