From 4fc05ea3382943c9a7e6ade6d0d7659209996237 Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Mon, 21 Jul 2025 11:19:27 -0400 Subject: [PATCH] new file: GainMatchSX3Front.C using a 2 step fit for the same method, going to implement a inverse fit now wherein the fronts are fit first and then the sum is fit against the backs. new file: GainMatchSX3Front.h new file: GainMatchSX3Front1.C --- .vscode/settings.json | 4 +- GainMatchSX3.C | 226 +++++++++++++--------------- GainMatchSX3Front.C | 338 ++++++++++++++++++++++++++++++++++++++++++ GainMatchSX3Front.h | 104 +++++++++++++ GainMatchSX3Front1.C | 245 ++++++++++++++++++++++++++++++ 5 files changed, 793 insertions(+), 124 deletions(-) create mode 100644 GainMatchSX3Front.C create mode 100644 GainMatchSX3Front.h create mode 100644 GainMatchSX3Front1.C diff --git a/.vscode/settings.json b/.vscode/settings.json index d94e58d..bbc9126 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -111,7 +111,9 @@ "RelBack_Fix_new.C": "cpp", "SiRelativeGains_Step1_new.C": "cpp", "charconv": "cpp", - "format": "cpp" + "format": "cpp", + "GainMatchSX3Front.C": "cpp", + "GainMatchSX3Front1.C": "cpp" }, "github-enterprise.uri": "https://fsunuc.physics.fsu.edu" } \ No newline at end of file diff --git a/GainMatchSX3.C b/GainMatchSX3.C index b98a28d..27b63f1 100644 --- a/GainMatchSX3.C +++ b/GainMatchSX3.C @@ -12,8 +12,6 @@ #include #include #include "Armory/ClassSX3.h" -#include "TGraphErrors.h" -#include "TMultiDimFit.h" #include "TVector3.h" @@ -28,10 +26,8 @@ int padID = 0; SX3 sx3_contr; TCutG *cut; -TCutG *cut1; std::map, std::vector>> dataPoints; - void GainMatchSX3::Begin(TTree * /*tree*/) { TString option = GetOption(); @@ -60,21 +56,6 @@ 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"); - if (!cutFile1 || cutFile1->IsZombie()) - { - std::cerr << "Error: Could not open UvD.root" << std::endl; - return; - } - 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"); } Bool_t GainMatchSX3::Process(Long64_t entry) @@ -104,11 +85,11 @@ Bool_t GainMatchSX3::Process(Long64_t entry) for (int i = 0; i < sx3.multi; i++) { - for (int j = i + 1; j < sx3.multi; j++) - { - if (sx3.id[i] == 3) - hsx3Coin->Fill(sx3.index[i], sx3.index[j]); - } + // for (int j = i + 1; j < sx3.multi; j++) + // { + // if (sx3.id[i] == 3) + // hsx3Coin->Fill(sx3.index[i], sx3.index[j]); + // } if (sx3.e[i] > 100) { ID.push_back(std::pair(sx3.id[i], i)); @@ -186,7 +167,7 @@ Bool_t GainMatchSX3::Process(Long64_t entry) for (int i = 0; i < sx3.multi; i++) { - if (sx3.id[i] == 3 && sx3.e[i] > 100) + if (sx3.id[i] == 3) { // Fill the histogram for the front vs back with gain correction hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk); @@ -208,7 +189,7 @@ Bool_t GainMatchSX3::Process(Long64_t entry) hist2d->Fill(sx3EUp + sx3EDn, sx3EBk); - if (cut && cut->IsInside(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)) { // Accumulate data for gain matching @@ -224,9 +205,6 @@ Bool_t GainMatchSX3::Process(Long64_t entry) void GainMatchSX3::Terminate() { - - // --- Store fit coefficients in memory --- - std::map, TVectorD> fitCoefficients; const int MAX_DET = 24; const int MAX_UP = 4; const int MAX_DOWN = 4; @@ -234,124 +212,126 @@ void GainMatchSX3::Terminate() double gainArray[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; bool gainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; + double fbgain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; + bool fbgainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; - std::ofstream outFile("sx3_MultiDimFit_results.txt"); + // std::map updn2DHistos; + std::map upCorrFactor; + + // === Gain matching === + + std::ofstream outFile("sx3_GainMatchback.txt"); if (!outFile.is_open()) { std::cerr << "Error opening output file!" << std::endl; return; } - // === Loop over all (id, bk, up, dn) combinations === - for (const auto &kv : dataPoints) { + // Gain fit using up+dn vs bk + for (const auto &kv : dataPoints) + + { + // kv.first is a tuple of (id, up, bk) + // kv.second is a vector of tuples (bkE, upE, dnE) auto [id, bk, u, d] = kv.first; const auto &pts = kv.second; + // Check if we have enough points for fitting + if (pts.size() < 5) + continue; - if (pts.size() < 20) continue; + std::vector bkE, udE; - std::vector x_bk, x_up, y_fsum; - for (const auto &pr : pts) { - double eBk, eUp, eDn; + for (const auto &pr : pts) + { + double eUp, eDn, eBk; std::tie(eBk, eUp, eDn) = pr; - if (eBk > 0 && eUp > 0 && eDn > 0) { - x_bk.push_back(eBk); - x_up.push_back(eUp); - y_fsum.push_back(eUp + eDn); + bkE.push_back(eBk); + udE.push_back(eUp + eDn); + } + + // Fill the TGraph with bkE and udE + 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) + { + continue; // Skip this fit if the slope is too far from 1 + } + 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) + { + for (int u = 0; u < MAX_UP; ++u) + { + 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; + } + } + } } } - - int nPoints = y_fsum.size(); - if (nPoints < 20) continue; - - TMultiDimFit *mdf = new TMultiDimFit(2, TMultiDimFit::kMonomials); - mdf->SetMaxPowers(new Int_t[2]{1, 1}); - mdf->SetMinAngle(10); - mdf->SetMinRelativeError(1e-4); - - double *x_row = new double[2]; - for (int i = 0; i < nPoints; ++i) { - x_row[0] = x_bk[i]; - x_row[1] = x_up[i]; - mdf->AddRow(x_row, y_fsum[i]); - } - delete[] x_row; - - mdf->Fit(); - - const TVectorD *coeffs = mdf->GetCoefficients(); - if (!coeffs || coeffs->GetNoElements() == 0 || !TMath::Finite((*coeffs)(0))) { - std::cerr << "Fit failed for Det" << id << " B" << bk << " U" << u << " D" << d << std::endl; - delete mdf; - continue; - } - - // Store coefficients in the map and write to file - fitCoefficients[kv.first] = *coeffs; - - int nCoeffs = mdf->GetNCoefficients(); - outFile << id << " " << bk << " " << u << " " << d; - printf("Fit for Det%d B%d U%d D%d -> ", id, bk, u, d); - for (int i = 0; i < nCoeffs; ++i) { - outFile << " " << (*coeffs)(i); - printf("p%d: %.4f ", i, (*coeffs)(i)); - } - outFile << std::endl; - printf("\n"); - - delete mdf; } + // for (int bk = 0; bk < MAX_BK; ++bk) + // { + // TString name = Form("hUpDnVsBk_%d", bk); + // TString title = Form("Up/Bk vs Dn/Bk for Back %d;Dn/Bk;Up/Bk", bk); + // updn2DHistos[bk] = new TH2F(name, title, 400, 0, 1, 400, 0, 1); + // } + outFile.close(); - std::cout << "Multi-dimensional gain matching complete. Results saved." << std::endl; + std::cout << "Gain matching complete." << std::endl; - // --- Stage 2: Apply corrections and create new histograms --- - std::cout << "--- Stage 2: Applying Corrections and Visualizing Results ---" << std::endl; - TH2F *hCorrectedFvB = new TH2F("hCorrectedFvB", "Gain Corrected Data;Predicted Front Sum (from fit);Measured Front Sum", 400, 0, 16000, 400, 0, 16000); + // === Create histograms === + TH2F *hFVB = new TH2F("hFVB", "Corrected Up+Dn vs Corrected Back;Corrected Back E;Up+Dn E", + 400, 0, 16000, 400, 0, 16000); + TH2F *hAsym = new TH2F("hAsym", "Up vs Dn dvide corrected back;Up/Back E;Dn/Back E", + 400, 0.0, 1.0, 400, 0.0, 1.0); - for (const auto &kv : dataPoints) { - // Find the coefficients for this segment - if (fitCoefficients.find(kv.first) == fitCoefficients.end()) { - continue; // Skip if no valid fit was found - } - const TVectorD& coeffs = fitCoefficients[kv.first]; - double p0 = coeffs(0); - double p1 = coeffs(1); - double p2 = coeffs(2); + // Fill histograms + for (const auto &kv : dataPoints) + { + auto [id, u, d, bk] = kv.first; + if (!gainValid[id][u][d][bk]) + continue; + double gain = gainArray[id][u][d][bk]; - // Loop over the data points for this segment and apply the correction - const auto &pts = kv.second; - for (const auto &pr : pts) { + // Prepare vectors to hold the points for TGraph + std::vector xVals; + std::vector yVals; + + for (const auto &pr : kv.second) + { double eBk, eUp, eDn; std::tie(eBk, eUp, eDn) = pr; - // Calculate the predicted front sum using the fit parameters - double predicted_front_sum = p0 + p1 * eBk + p2 * eUp; - - // The measured front sum is just the raw sum - double measured_front_sum = eUp + eDn; + double updn = eUp + eDn; + if (updn == 0 || eBk == 0) + continue; - // Fill the corrected histogram - hCorrectedFvB->Fill(predicted_front_sum, measured_front_sum); + double asym = (eUp - eDn) / updn; + double correctedBack = eBk * gain; + + hFVB->Fill(correctedBack, updn); + hAsym->Fill(eUp / correctedBack, eDn / correctedBack); } } - - // --- Stage 3: Draw the comparison canvases --- - 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(2); - hCorrectedFvB->SetTitle("After Correction"); - hCorrectedFvB->Draw("colz"); - // Draw a perfect y=x line for comparison - TF1 *diag = new TF1("diag", "x", 0, 16000); - diag->SetLineColor(kRed); - diag->SetLineWidth(2); - diag->Draw("same"); } \ No newline at end of file diff --git a/GainMatchSX3Front.C b/GainMatchSX3Front.C new file mode 100644 index 0000000..142dcd5 --- /dev/null +++ b/GainMatchSX3Front.C @@ -0,0 +1,338 @@ +#define GainMatchSX3Front_cxx + +#include "GainMatchSX3Front.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "Armory/ClassSX3.h" +#include "TGraphErrors.h" +#include "TMultiDimFit.h" + +#include "TVector3.h" + +TH2F *hSX3FvsB; +TH2F *hSX3FvsB_g; +TH2F *hsx3IndexVE; +TH2F *hsx3IndexVE_g; +TH2F *hSX3; +TH2F *hsx3Coin; + +int padID = 0; + +SX3 sx3_contr; +TCutG *cut; +TCutG *cut1; +std::map, std::vector>> dataPoints; + +// Gain arrays + +const int MAX_DET = 24; +const int MAX_UP = 4; +const int MAX_DOWN = 4; +const int MAX_BK = 4; +double backGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; +bool backGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; +double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; +bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; + +void GainMatchSX3Front::Begin(TTree * /*tree*/) +{ + TString option = GetOption(); + + hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000); + hSX3FvsB_g = new TH2F("hSX3FvsB_g", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000); + hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000); + hsx3IndexVE_g = new TH2F("hsx3IndexVE_g", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000); + hSX3 = new TH2F("hSX3", "SX3 Front v Back; Fronts; Backs", 8, 0, 8, 4, 0, 4); + + hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12); + + sx3_contr.ConstructGeo(); + + // Load the TCutG object + TFile *cutFile = TFile::Open("sx3cut.root"); + bool cutLoaded = (cut != nullptr); + cut = dynamic_cast(cutFile->Get("sx3cut")); + if (!cut) + { + std::cerr << "Error: Could not find TCutG named 'sx3cut' in sx3cut.root" << std::endl; + 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_GainMatchback.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) + { + backGain[id][bk][u][d] = gain; + backGainValid[id][bk][u][d] = true; + } + + infile.close(); + std::cout << "Loaded back gains from " << filename << std::endl; + SX3 sx3_contr; +} + +Bool_t GainMatchSX3Front::Process(Long64_t entry) +{ + + b_sx3Multi->GetEntry(entry); + b_sx3ID->GetEntry(entry); + b_sx3Ch->GetEntry(entry); + b_sx3E->GetEntry(entry); + b_sx3T->GetEntry(entry); + + sx3.CalIndex(); + + std::vector> ID; + for (int i = 0; i < sx3.multi; i++) + { + + for (int j = i + 1; j < sx3.multi; j++) + { + if (sx3.id[i] == 3) + hsx3Coin->Fill(sx3.index[i], sx3.index[j]); + } + if (sx3.e[i] > 100) + { + ID.push_back(std::pair(sx3.id[i], i)); + hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]); + } + } + + if (ID.size() > 0) + { + std::sort(ID.begin(), ID.end(), [](const std::pair &a, const std::pair &b) + { return a.first < b.first; }); + + // start with the first entry in the sorted array: channels that belong to the same detector are together in sequenmce + std::vector> sx3ID; + sx3ID.push_back(ID[0]); + bool found = false; + + for (size_t i = 1; i < ID.size(); i++) + { // Check if id of i belongs to the same detector and then add it to the detector ID vector + if (ID[i].first == sx3ID.back().first) + { // count the nunmber of hits that belong to the same detector + sx3ID.push_back(ID[i]); + + if (sx3ID.size() >= 3) + { + found = true; + } + } + else + { // the next event does not belong to the same detector, abandon the first event and continue with the next one + if (!found) + { + sx3ID.clear(); + sx3ID.push_back(ID[i]); + } + } + } + + if (found) + { + int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1; + float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0; + + 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.ch[index] % 2 == 0) + { + sx3ChDn = sx3.ch[index] / 2; + sx3EDn = sx3.e[index]; + } + else + { + sx3ChUp = sx3.ch[index] / 2; + sx3EUp = 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]; + } + } + // If we have a valid front and back channel, fill the histograms + hSX3->Fill(sx3ChDn + 4, sx3ChBk); + hSX3->Fill(sx3ChUp, sx3ChBk); + + // Fill the histogram for the front vs back + hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk); + + for (int i = 0; i < sx3.multi; i++) + { + if (sx3.id[i] == 3 && sx3.e[i] > 100) + { + // back gain correction + + // Fill the histogram for the front vs back with gain correction + hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk); + // Fill the index vs energy histogram + hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]); + // } + // { + TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk); + TH2F *hist2d = (TH2F *)gDirectory->Get(histName); + if (!hist2d) + { + hist2d = new TH2F(histName, Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk), 400, 0, 16000, 400, 0, 16000); + } + + // if (sx3ChBk == 2) + // printf("Found back channel Det %d Back %d \n", sx3.id[i], sx3ChBk); + // hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]); + // hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk); + + hist2d->Fill(sx3EUp + sx3EDn, sx3EBk); + + if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk) && + cut1 && cut1->IsInside(sx3EUp / sx3EBk, sx3EDn / sx3EBk)) + { + + if (backGainValid[sx3.id[i]][sx3ChBk][sx3ChUp][sx3ChDn]) + { + sx3EBk *= backGain[sx3.id[i]][sx3ChBk][sx3ChUp][sx3ChDn]; + } + // 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); + } + } + } + } + } + + return kTRUE; +} + +void GainMatchSX3Front::Terminate() +{ + + std::map, TVectorD> fitCoefficients; + + // === Gain matching === + + std::ofstream outFile("sx3_GainMatchfront.txt"); + if (!outFile.is_open()) + { + std::cerr << "Error opening output file!" << std::endl; + return; + } + + TH2F *hUvD = new TH2F("hUvD", " UvD; Up/CorrBack; Down/CorrBack", 600, 0, 1, 600, 0, 1); + + for (const auto &kv : dataPoints) + { + auto [id, bk, u, d] = kv.first; + const auto &pts = kv.second; + + if (pts.size() < 5) + continue; + + std::vector udE, corrBkE; + + for (const auto &pr : pts) + { + double eBkCorr, eUp, eDn; + std::tie(eBkCorr, eUp, eDn) = pr; + udE.push_back(eUp + eDn); + corrBkE.push_back(eBkCorr); + hUvD->Fill(eUp / eBkCorr, eDn / eBkCorr); + } + + 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; + + outFile << id << " " << bk << " " << u << " " << d << " " << frontGain[id][bk][u][d] << std::endl; + printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]); + } + + outFile.close(); + 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 *hCorrectedUvD = new TH2F("hCorrectedUvD", "Corrected UvD; UvD Up; UvD Down", 600, 0, 1, 600, 0, 1); + + for (const auto &kv : dataPoints) + { + + auto [id, bk, u, d] = kv.first; + double front = frontGain[id][bk][u][d]; + + for (const auto &pr : kv.second) + { + double eBk, eUp, eDn; + std::tie(eBk, eUp, eDn) = pr; + double corrUp = eUp * front; + double corrDn = eDn * front; + + hCorrectedFvB->Fill(corrUp + corrDn, eBk); + hCorrectedUvD->Fill(corrUp / eBk, corrDn / eBk); + } + } + + // === 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(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 diff --git a/GainMatchSX3Front.h b/GainMatchSX3Front.h new file mode 100644 index 0000000..47aef61 --- /dev/null +++ b/GainMatchSX3Front.h @@ -0,0 +1,104 @@ +#ifndef GainMatchSX3Front_h +#define GainMatchSX3Front_h + +#include +#include +#include +#include + +#include "Armory/ClassDet.h" + +class GainMatchSX3Front : public TSelector { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + +// Fixed size dimensions of array or collections stored in the TTree if any. + + // Declaration of leaf types + Det sx3; + Det qqq; + Det pc ; + + ULong64_t evID; + UInt_t run; + + // List of branches + TBranch *b_eventID; //! + TBranch *b_run; //! + TBranch *b_sx3Multi; //! + TBranch *b_sx3ID; //! + TBranch *b_sx3Ch; //! + TBranch *b_sx3E; //! + TBranch *b_sx3T; //! + TBranch *b_qqqMulti; //! + TBranch *b_qqqID; //! + TBranch *b_qqqCh; //! + TBranch *b_qqqE; //! + TBranch *b_qqqT; //! + TBranch *b_pcMulti; //! + TBranch *b_pcID; //! + TBranch *b_pcCh; //! + TBranch *b_pcE; //! + TBranch *b_pcT; //! + + GainMatchSX3Front(TTree * /*tree*/ =0) : fChain(0) { } + virtual ~GainMatchSX3Front() { } + virtual Int_t Version() const { return 2; } + virtual void Begin(TTree *tree); + virtual void SlaveBegin(TTree *tree); + virtual void Init(TTree *tree); + virtual Bool_t Notify(); + virtual Bool_t Process(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; } + virtual void SetOption(const char *option) { fOption = option; } + virtual void SetObject(TObject *obj) { fObject = obj; } + virtual void SetInputList(TList *input) { fInput = input; } + virtual TList *GetOutputList() const { return fOutput; } + virtual void SlaveTerminate(); + virtual void Terminate(); + + ClassDef(GainMatchSX3Front,0); +}; + +#endif + +#ifdef GainMatchSX3Front_cxx +void GainMatchSX3Front::Init(TTree *tree){ + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("evID", &evID, &b_eventID); + fChain->SetBranchAddress("run", &run, &b_run); + + sx3.SetDetDimension(24,12); + qqq.SetDetDimension(4,32); + pc.SetDetDimension(2,24); + + fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi); + fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID); + fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch); + fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E); + fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T); + +} + +Bool_t GainMatchSX3Front::Notify(){ + + return kTRUE; +} + +void GainMatchSX3Front::SlaveBegin(TTree * /*tree*/){ + + TString option = GetOption(); + +} + +void GainMatchSX3Front::SlaveTerminate(){ + +} + + +#endif // #ifdef GainMatchSX3Front_cxx \ No newline at end of file diff --git a/GainMatchSX3Front1.C b/GainMatchSX3Front1.C new file mode 100644 index 0000000..28ac0e2 --- /dev/null +++ b/GainMatchSX3Front1.C @@ -0,0 +1,245 @@ +#define GainMatchSX3_cxx + +#include "GainMatchSX3.h" +#include "Armory/ClassSX3.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Constants +const int MAX_DET = 24; +const int MAX_BK = 4; +const int MAX_UP = 4; +const int MAX_DOWN = 4; + +// Gain arrays +double backGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; +bool backGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; + +double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; +bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; + +// Data container +std::map, std::vector>> dataPoints; + +// Load back gains +void LoadBackGains(const std::string &filename) +{ + 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) + { + backGain[id][bk][u][d] = gain; + backGainValid[id][bk][u][d] = true; + } + + infile.close(); + std::cout << "Loaded back gains from " << filename << std::endl; + SX3 sx3_contr; +} + +// Front gain matching function +Bool_t GainMatchSX3::Process(Long64_t entry) +{ + // Link SX3 branches + + b_sx3Multi->GetEntry(entry); + b_sx3ID->GetEntry(entry); + b_sx3Ch->GetEntry(entry); + b_sx3E->GetEntry(entry); + b_sx3T->GetEntry(entry); + + sx3.CalIndex(); + + Long64_t nentries = tree->GetEntries(Long64_t entry); + std::cout << "Total entries: " << nentries << std::endl; + + TH2F *hBefore = new TH2F("hBefore", "Before Correction;E_Up+E_Dn;Back Energy", 400, 0, 40000, 400, 0, 40000); + TH2F *hAfter = new TH2F("hAfter", "After Correction;E_Up+E_Dn;Corrected Back Energy", 400, 0, 40000, 400, 0, 40000); + + for (Long64_t entry = 0; entry < nentries; ++entry) + { + tree->GetEntry(entry); + sx3.CalIndex(); + + std::vector> ID; + + for (int i = 0; i < sx3.multi; i++) + { + if (sx3.e[i] > 100) + { + ID.push_back({sx3.id[i], i}); + } + } + + if (ID.empty()) + continue; + + // Sort by id + std::sort(ID.begin(), ID.end(), [](auto &a, auto &b) { return a.first < b.first; }); + + std::vector> sx3ID; + sx3ID.push_back(ID[0]); + bool found = false; + + for (size_t i = 1; i < ID.size(); i++) + { + if (ID[i].first == sx3ID.back().first) + { + sx3ID.push_back(ID[i]); + if (sx3ID.size() >= 3) + found = true; + } + else if (!found) + { + sx3ID.clear(); + sx3ID.push_back(ID[i]); + } + } + + if (!found) + continue; + + int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1; + float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0; + int sx3id = sx3ID[0].first; + + for (auto &[id, idx] : sx3ID) + { + if (sx3.ch[idx] < 8) + { + if (sx3.ch[idx] % 2 == 0) + { + sx3ChDn = sx3.ch[idx] / 2; + sx3EDn = sx3.e[idx]; + } + else + { + sx3ChUp = sx3.ch[idx] / 2; + sx3EUp = sx3.e[idx]; + } + } + else + { + sx3ChBk = sx3.ch[idx] - 8; + sx3EBk = sx3.e[idx]; + } + } + + if (sx3ChUp < 0 || sx3ChDn < 0 || sx3ChBk < 0) + continue; + + if (!backGainValid[sx3id][sx3ChBk][sx3ChUp][sx3ChDn]) + continue; + + double corrBk = sx3EBk * backGain[sx3id][sx3ChBk][sx3ChUp][sx3ChDn]; + + hBefore->Fill(sx3EUp + sx3EDn, sx3EBk); + hAfter->Fill(sx3EUp + sx3EDn, corrBk); + + dataPoints[{sx3id, sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(corrBk, sx3EUp, sx3EDn); + } + + // === Fit front gains === + std::ofstream outFile("sx3_GainMatchfront.txt"); + if (!outFile.is_open()) + { + std::cerr << "Error opening sx3_GainMatchfront.txt!" << std::endl; + return; + } + + for (const auto &kv : dataPoints) + { + auto [id, bk, u, d] = kv.first; + const auto &pts = kv.second; + + if (pts.size() < 5) + continue; + + std::vector udE, corrBkE; + + for (const auto &pr : pts) + { + double eBkCorr, eUp, eDn; + std::tie(eBkCorr, eUp, eDn) = pr; + udE.push_back(eUp + eDn); + corrBkE.push_back(eBkCorr); + } + + 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; + + outFile << id << " " << bk << " " << u << " " << d << " " << frontGain[id][bk][u][d] << std::endl; + printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]); + } + + outFile.close(); + std::cout << "Front gain matching complete." << std::endl; + + // === Draw diagnostic plots === + gStyle->SetOptStat(1110); + TCanvas *c = new TCanvas("c", "Gain Matching Diagnostics", 1200, 600); + c->Divide(2, 1); + + c->cd(1); + hBefore->Draw("colz"); + TF1 *diag1 = new TF1("diag1", "x", 0, 40000); + diag1->SetLineColor(kRed); + diag1->Draw("same"); + + c->cd(2); + hAfter->Draw("colz"); + TF1 *diag2 = new TF1("diag2", "x", 0, 40000); + diag2->SetLineColor(kRed); + diag2->Draw("same"); +} + +int main(int argc, char **argv) +{ + TApplication app("app", &argc, argv); + + // Load back gains + LoadBackGains("sx3_GainMatchback.txt"); + + // Open tree + TFile *f = TFile::Open("input_tree.root"); // <<< Change file name + if (!f || f->IsZombie()) + { + std::cerr << "Cannot open input_tree.root!" << std::endl; + return 1; + } + TTree *tree = (TTree *)f->Get("tree"); + if (!tree) + { + std::cerr << "Tree not found!" << std::endl; + return 1; + } + + // Run front gain matching + GainMatchSX3(tree); + + app.Run(); + return 0; +}