From 3d0d176f5a7392e3d125d594788346049a315ad9 Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Thu, 24 Jul 2025 16:12:46 -0400 Subject: [PATCH] modified: GainMatchSX3.C 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. --- GainMatchSX3.C | 189 +++++++++++++++++++++++++++++++++----------- GainMatchSX3Front.C | 117 +++++++++++++++++++-------- 2 files changed, 225 insertions(+), 81 deletions(-) diff --git a/GainMatchSX3.C b/GainMatchSX3.C index 27b63f1..366b176 100644 --- a/GainMatchSX3.C +++ b/GainMatchSX3.C @@ -12,7 +12,7 @@ #include #include #include "Armory/ClassSX3.h" - +#include #include "TVector3.h" TH2F *hSX3FvsB; @@ -26,7 +26,17 @@ int padID = 0; SX3 sx3_contr; TCutG *cut; +TCutG *cut1; std::map, std::vector>> dataPoints; +std::map, 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(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 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; } } } diff --git a/GainMatchSX3Front.C b/GainMatchSX3Front.C index 142dcd5..941380d 100644 --- a/GainMatchSX3Front.C +++ b/GainMatchSX3Front.C @@ -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 udE, corrBkE; + std::vector 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 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; } \ No newline at end of file