From 265ebd3372d845d5f308f2446b89cb0c3b7b979b Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Tue, 30 Sep 2025 15:17:00 -0400 Subject: [PATCH] modified: Calibration.C modified: GainMatchQQQ.C modified: GainMatchSX3.C modified: GainMatchSX3Front.C trying out not doing back gainmatching to see if that improves fits. --- Calibration.C | 351 +++++++++++++++++++++++++------------------- GainMatchQQQ.C | 86 +---------- GainMatchSX3.C | 202 ++++++++++--------------- GainMatchSX3Front.C | 8 +- 4 files changed, 284 insertions(+), 363 deletions(-) diff --git a/Calibration.C b/Calibration.C index a0c3a41..68e87f7 100644 --- a/Calibration.C +++ b/Calibration.C @@ -53,13 +53,15 @@ const int MAX_BK = 4; const int MAX_QQQ = 4; const int MAX_RING = 16; const int MAX_WEDGE = 16; -double backGain[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; -bool backGainValid[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; +// double backGain[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; +// bool backGainValid[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; double frontGain[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; bool frontGainValid[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; double uvdslope[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; double qqqGain[MAX_QQQ][MAX_BK][MAX_UP] = {{{0}}}; bool qqqGainValid[MAX_QQQ][MAX_BK][MAX_UP] = {{{false}}}; +TH1F *hSX3Spectra[MAX_SX3][MAX_BK][MAX_UP][MAX_DOWN]; +TH1F *hQQQSpectra[MAX_QQQ][MAX_RING][MAX_WEDGE]; void Calibration::Begin(TTree * /*tree*/) { @@ -82,26 +84,26 @@ void Calibration::Begin(TTree * /*tree*/) sx3_contr.ConstructGeo(); pw_contr.ConstructGeo(); // ----------------------- Load Back Gains - { - std::string filename = "sx3_GainMatchback.txt"; - std::ifstream infile(filename); - if (!infile.is_open()) - { - std::cerr << "Error opening " << filename << "!" << std::endl; - } - else - { - 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] = (gain > 0); - } - infile.close(); - std::cout << "Loaded back gains from " << filename << std::endl; - } - } + // { + // std::string filename = "sx3_GainMatchback.txt"; + // std::ifstream infile(filename); + // if (!infile.is_open()) + // { + // std::cerr << "Error opening " << filename << "!" << std::endl; + // } + // else + // { + // 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] = (gain > 0); + // } + // infile.close(); + // std::cout << "Loaded back gains from " << filename << std::endl; + // } + // } // ----------------------- Load Front Gains { @@ -147,19 +149,42 @@ void Calibration::Begin(TTree * /*tree*/) } } + for (int id = 0; id < MAX_SX3; id++) + { + for (int bk = 0; bk < MAX_BK; bk++) + { + for (int up = 0; up < MAX_UP; up++) + { + for (int dn = 0; dn < MAX_DOWN; dn++) + { + TString hname = Form("hCal_id%d_bk%d_up%d_dn%d", id, bk, up, dn); + TString htitle = Form("SX3 id%d bk%d up%d dn%d; Energy (arb); Counts", id, bk, up, dn); + hSX3Spectra[id][bk][up][dn] = new TH1F(hname, htitle, 4000, 0, 16000); + } + } + } + } + for (int det = 0; det < MAX_QQQ; det++) + { + for (int ring = 0; ring < MAX_RING; ring++) + { + for (int wedge = 0; wedge < MAX_WEDGE; wedge++) + { + TString hname = Form("hCal_qqq%d_ring%d_wedge%d", det, ring, wedge); + TString htitle = Form("QQQ det%d ring%d wedge%d; Energy (arb); Counts", det, ring, wedge); + hQQQSpectra[det][ring][wedge] = new TH1F(hname, htitle, 4000, 0, 16000); + } + } + } + SX3 sx3_contr; } Bool_t Calibration::Process(Long64_t entry) { - - // if ( entry > 100 ) return kTRUE; - hitPos.Clear(); HitNonZero = false; - // if( entry > 1) return kTRUE; - // printf("################### ev : %llu \n", entry); - + // Load branches b_sx3Multi->GetEntry(entry); b_sx3ID->GetEntry(entry); b_sx3Ch->GetEntry(entry); @@ -177,46 +202,38 @@ Bool_t Calibration::Process(Long64_t entry) qqq.CalIndex(); pc.CalIndex(); - // sx3.Print(); - // ########################################################### Raw data - // //======================= SX3 sx3ecut = false; std::vector> ID; // first = id, 2nd = index for (int i = 0; i < sx3.multi; i++) { - ID.push_back(std::pair(sx3.id[i], i)); + ID.emplace_back(sx3.id[i], i); hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]); if (sx3.e[i] > 100) - { sx3ecut = true; - } for (int j = i + 1; j < sx3.multi; j++) - { hsx3Coin->Fill(sx3.index[i], sx3.index[j]); - } } - // --- safe SX3 handling (replace your existing block that builds sx3ID) --- - if (ID.size() > 0) + // --- SX3 safe handling --- + if (!ID.empty()) { - std::sort(ID.begin(), ID.end(), [](const std::pair &a, const std::pair &b) + 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 { @@ -224,141 +241,101 @@ Bool_t Calibration::Process(Long64_t entry) { sx3ID.clear(); sx3ID.push_back(ID[i]); - found = false; } } } - if (found) { - // initialize to sentinel values - int sx3ChUp = -1; - int sx3ChDn = -1; - int sx3ChBk = -1; - float sx3EUp = 0.0f; - float sx3EDn = 0.0f; - float sx3EBk = 0.0f; + int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1; + float sx3EUp = 0.0f, sx3EDn = 0.0f, sx3EBk = 0.0f; - // collect channels/energies - for (size_t i = 0; i < sx3ID.size(); i++) + for (auto &p : sx3ID) { - int index = sx3ID[i].second; + int index = p.second; int ch = sx3.ch[index]; float e = sx3.e[index]; - if (ch < 8) // front channels + if (ch < 8) { - // you used even/odd to denote down/up — keep that convention - if ((ch % 2) == 0) // down + if ((ch % 2) == 0) // even -> down { sx3ChDn = ch; sx3EDn = e; } - else // up + else // odd -> up { sx3ChUp = ch; sx3EUp = e; } } - else // back channels (assuming back channels are 8..11 or so) + else { - sx3ChBk = ch; // store as raw channel number; adapt if you index bk differently - sx3EBk = e; // if you want to track back energy too + sx3ChBk = ch; + sx3EBk = e; } } - // Basic sanity checks before using indices: - bool haveFrontPair = (sx3ChUp >= 0 && sx3ChDn >= 0); + bool haveFrontPair = (sx3ChUp >= 0 || sx3ChDn >= 0); bool haveBack = (sx3ChBk >= 0); + int sx3Id = sx3ID[0].first; - // convert raw channel numbers to array indices if needed: - int bk_index = (haveBack ? (sx3ChBk - 8) : -1); - int up_index = (haveFrontPair ? sx3ChUp : -1); - int dn_index = (haveFrontPair ? sx3ChDn : -1); - auto sx3Id = sx3ID[0].first; + // CORRECTED: map channel (0..7) to front-index (0..3) + int bk_index = (haveBack ? sx3ChBk - 8 : -1); + int up_index = (sx3ChUp >= 0 ? sx3ChUp / 2 : -1); // <<-- IMPORTANT FIX + int dn_index = (sx3ChDn >= 0 ? sx3ChDn / 2 : -1); // <<-- IMPORTANT FIX - double calibEUp, calibEDn, calibEBack = 0.0; + double GM_EUp = 0.0, GM_EDn = 0.0, calibEBack = 0.0; - if (haveFrontPair && haveBack) + if (haveBack) { - // If you stored front gains indexed by [id][bk][up][down] + // --- ALWAYS fill raw ADC for diagnostics + // (temporarily use the existing spectrum to confirm fills) + // If you don't want raw values mixed with calibrated later, create a separate _raw array. + hSX3Spectra[sx3Id][bk_index][up_index][dn_index]->Fill(sx3EUp); + + // --- If gain is available, also fill calibrated energy if (frontGainValid[sx3Id][bk_index][up_index][dn_index]) { - calibEUp = frontGain[sx3Id][bk_index][up_index][dn_index] * sx3EUp; - // calibEDn = frontGain[sx3Id][bk_index][up_index][dn_index] * sx3EDn; - } - if (backGainValid[sx3Id][bk_index][up_index][dn_index]) - { - calibEBack = backGain[sx3Id][bk_index][up_index][dn_index] * sx3EBk; - } - } - - // Only call CalSX3Pos if we have reasonable energies (avoid calling with zeros/uninitialized) - if (haveFrontPair && (calibEUp > 50.0) && haveBack && (calibEBack > 50.0)) - { - // find exact back energy value from sx3 entries if you tracked it above - float backEnergyRaw = 0.0f; - // locate the back index in sx3ID if needed - for (size_t k = 0; k < sx3ID.size(); ++k) - { - int idx = sx3ID[k].second; - if (sx3.ch[idx] >= 8) - { - backEnergyRaw = sx3.e[idx]; - break; - } + GM_EUp = frontGain[sx3Id][bk_index][up_index][dn_index] * sx3EUp; + if (GM_EUp > 50.0) + hSX3Spectra[sx3Id][bk_index][up_index][dn_index]->Fill(GM_EUp); // optional: mixes raw+calib } - hsx3IndexVE_gm->Fill(sx3.index[sx3ID[0].second], calibEUp); + // Keep the other diagnostic plots + hsx3IndexVE_gm->Fill(sx3.index[sx3ID[0].second], GM_EUp); hSX3->Fill(sx3ChDn + 4, sx3ChBk); hSX3->Fill(sx3ChUp, sx3ChBk); + hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk); - // Fill the histogram for the front vs back - hSX3FvsB->Fill(sx3EUp + sx3EDn, calibEBack); - - sx3_contr.CalSX3Pos(sx3Id, sx3ChUp, sx3ChDn, sx3ChBk, static_cast(calibEUp), static_cast(calibEDn)); - hitPos = sx3_contr.GetHitPos(); - HitNonZero = true; + if (GM_EUp > 50.0 && sx3EBk > 50.0) + { + sx3_contr.CalSX3Pos(sx3Id, sx3ChUp, sx3ChDn, sx3ChBk, GM_EUp, sx3EDn); + hitPos = sx3_contr.GetHitPos(); + HitNonZero = true; + } } - } // found + else + { + // Debug print for channels that didn't pass validation -- helps find indexing problems + static int dbgCount = 0; + if (dbgCount < 20) // only print first few to avoid flood + { + std::cout << Form("DEBUG SX3 skip: id=%d chUp=%d chDn=%d chBk=%d -> up_idx=%d dn_idx=%d bk_idx=%d", + sx3Id, sx3ChUp, sx3ChDn, sx3ChBk, up_index, dn_index, bk_index) + << std::endl; + dbgCount++; + } + } + } } - // //======================= QQQ + // ======================= QQQ ======================= for (int i = 0; i < qqq.multi; i++) { - - int det = qqq.id[i]; // detector ID (0–3) - int ch = qqq.ch[i]; // raw channel (0–31) - - // Separate ring vs wedge channel - int ring = -1; - int wedge = -1; - if (ch < 16) - { // wedge - wedge = ch; - } - else - { // ring - ring = ch - 16; - } - - double Ecal = qqq.e[i]; // default = raw - if (ring >= 0 && wedge >= 0 && qqqGainValid[det][ring][wedge]) - { - Ecal *= qqqGain[det][ring][wedge]; - } - // for( int j = 0; j < pc.multi; j++){ - // if(pc.index[j]==4){ - hqqqIndexVE_gm->Fill(qqq.index[i], Ecal); - hqqqIndexVE->Fill(qqq.index[i], qqq.e[i]); - - // } - // printf("QQQ ID : %d, ch : %d, e : %d \n", qqq.id[i], qqq.ch[i], qqq.e[i]); + int det = qqq.id[i]; if (qqq.e[i] > 100) - { qqqEcut = true; - } - // } + for (int j = 0; j < qqq.multi; j++) { if (j == i) @@ -366,18 +343,11 @@ Bool_t Calibration::Process(Long64_t entry) hqqqCoin->Fill(qqq.index[i], qqq.index[j]); } - // } - for (int j = i + 1; j < qqq.multi; j++) { - // if( qqq.used[i] == true ) continue; - - // if( qqq.id[i] == qqq.id[j] && (16 - qqq.ch[i]) * (16 - qqq.ch[j]) < 0 ){ // must be same detector and wedge and ring if (qqq.id[i] == qqq.id[j]) - { // must be same detector - - int chWedge = -1; - int chRing = -1; + { + int chWedge = -1, chRing = -1; if (qqq.ch[i] < qqq.ch[j]) { chRing = qqq.ch[j] - 16; @@ -388,15 +358,29 @@ Bool_t Calibration::Process(Long64_t entry) chRing = qqq.ch[i]; chWedge = qqq.ch[j] - 16; } - // printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing); + + double Ecal = qqq.e[i]; + if (det >= 0 && det < MAX_QQQ && + chRing >= 0 && chRing < MAX_RING && + chWedge >= 0 && chWedge < MAX_WEDGE) + { + // ALWAYS fill raw energy for diagnostics + hQQQSpectra[det][chRing][chWedge]->Fill(qqq.e[i]); + + // If calibrated gain is present, also fill calibrated energy + if (qqqGainValid[det][chRing][chWedge]) + { + double Ecal = qqq.e[i] * qqqGain[det][chRing][chWedge]; + hQQQSpectra[det][chRing][chWedge]->Fill(Ecal); // optional: mixes raw+calib + } + } + + hqqqIndexVE_gm->Fill(qqq.index[i], Ecal); + hqqqIndexVE->Fill(qqq.index[i], qqq.e[i]); double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5); double rho = 50. + 40. / 16. * (chRing + 0.5); - // if(qqq.e[i]>50){ hqqqPolar->Fill(theta, rho); - // } - // qqq.used[i] = true; - // qqq.used[j] = true; if (!HitNonZero) { @@ -413,4 +397,75 @@ Bool_t Calibration::Process(Long64_t entry) } void Calibration::Terminate() { -} \ No newline at end of file + const double AM241_ALPHA = 5486.0; // keV + + // ----------------------- Summary Plots + TH2F *hSX3Summary = new TH2F("hSX3Summary", "SX3 Channel Means;Channel Index;Mean (ADC)", + MAX_SX3 * MAX_BK * MAX_UP * MAX_DOWN, 0, MAX_SX3 * MAX_BK * MAX_UP * MAX_DOWN, + 200, 0, 10000); + + TH2F *hQQQSummary = new TH2F("hQQQSummary", "QQQ Channel Means;Channel Index;Mean (ADC)", + MAX_QQQ * MAX_RING * MAX_WEDGE, 0, MAX_QQQ * MAX_RING * MAX_WEDGE, + 200, 0, 10000); + + // ----------------------- SX3 Calibration (quick check with mean) + for (int id = 0; id < MAX_SX3; id++) + { + for (int bk = 0; bk < MAX_BK; bk++) + { + for (int up = 0; up < MAX_UP; up++) + { + for (int dn = 0; dn < MAX_DOWN; dn++) + { + TH1F *hSpec = hSX3Spectra[id][bk][up][dn]; + if (!hSpec || hSpec->GetEntries() < 200) + continue; + + double mean = hSpec->GetMean(); + + int sx3Index = (((id * MAX_BK + bk) * MAX_UP + up) * MAX_DOWN + dn); + hSX3Summary->Fill(sx3Index, mean); + + std::cout << Form("SX3 id%d bk%d up%d dn%d → mean %.1f", + id, bk, up, dn, mean) + << std::endl; + } + } + } + } + + // ----------------------- QQQ Calibration (quick check with mean) + for (int det = 0; det < MAX_QQQ; det++) + { + for (int ring = 0; ring < MAX_RING; ring++) + { + for (int wedge = 0; wedge < MAX_WEDGE; wedge++) + { + TH1F *hSpec = hQQQSpectra[det][ring][wedge]; + if (!hSpec || hSpec->GetEntries() < 200) + continue; + + double mean = hSpec->GetMean(); + + int qqqIndex = ((det * MAX_RING + ring) * MAX_WEDGE + wedge); + hQQQSummary->Fill(qqqIndex, mean); + + std::cout << Form("QQQ det%d ring%d wedge%d → mean %.1f", + det, ring, wedge, mean) + << std::endl; + } + } + } + + // ----------------------- Draw Summary + TCanvas *cSum = new TCanvas("cSum", "Calibration Summary (Means)", 1200, 600); + cSum->Divide(2, 1); + + cSum->cd(1); + hSX3Summary->Draw("COLZ"); + + cSum->cd(2); + hQQQSummary->Draw("COLZ"); + + cSum->Update(); +} diff --git a/GainMatchQQQ.C b/GainMatchQQQ.C index 814b25c..c640fe1 100644 --- a/GainMatchQQQ.C +++ b/GainMatchQQQ.C @@ -11,16 +11,13 @@ #include #include -#include "Armory/ClassSX3.h" #include "TVector3.h" -TH2F *hSX3FvsB; TH2F *hQQQFVB; int padID = 0; -SX3 sx3_contr; TCutG *cut; std::map, std::vector>> dataPoints; @@ -28,10 +25,8 @@ void GainMatchQQQ::Begin(TTree * /*tree*/) { TString option = GetOption(); - hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000); - hQQQFVB = new TH2F("hQQQFVB", "number of good QQQ vs QQQ id", 400, 0, 16000, 400, 0, 16000); + hQQQFVB = new TH2F("hQQQFVB", "QQQ Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000); - sx3_contr.ConstructGeo(); // Load the TCutG object TFile *cutFile = TFile::Open("qqqcorr.root"); @@ -51,92 +46,13 @@ void GainMatchQQQ::Begin(TTree * /*tree*/) Bool_t GainMatchQQQ::Process(Long64_t entry) { - - b_sx3Multi->GetEntry(entry); - b_sx3ID->GetEntry(entry); - b_sx3Ch->GetEntry(entry); - b_sx3E->GetEntry(entry); - b_sx3T->GetEntry(entry); b_qqqMulti->GetEntry(entry); b_qqqID->GetEntry(entry); b_qqqCh->GetEntry(entry); b_qqqE->GetEntry(entry); b_qqqT->GetEntry(entry); - b_pcMulti->GetEntry(entry); - b_pcID->GetEntry(entry); - b_pcCh->GetEntry(entry); - b_pcE->GetEntry(entry); - b_pcT->GetEntry(entry); - sx3.CalIndex(); qqq.CalIndex(); - pc.CalIndex(); - - std::vector> ID; - for (int i = 0; i < sx3.multi; i++) - { - ID.push_back(std::pair(sx3.id[i], i)); - } - - if (ID.size() > 0) - { - std::sort(ID.begin(), ID.end(), [](const std::pair &a, const std::pair &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) - { - 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; - if (sx3.ch[index] < 8) - { - if (sx3.ch[index] % 2 == 0) - { - sx3ChDn = sx3.ch[index]; - sx3EDn = sx3.e[index]; - } - else - { - sx3ChUp = sx3.ch[index]; - sx3EUp = sx3.e[index]; - } - } - else - { - sx3ChBk = sx3.ch[index]; - sx3EBk = sx3.e[index]; - } - } - hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk); - } - } for (int i = 0; i < qqq.multi; i++) { diff --git a/GainMatchSX3.C b/GainMatchSX3.C index 35a179a..483cbf3 100644 --- a/GainMatchSX3.C +++ b/GainMatchSX3.C @@ -35,14 +35,13 @@ 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}}}}; +double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; +bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; // ==== Configuration Flags ==== const bool interactiveMode = false; // If true: show canvas + wait for user const bool verboseFit = true; // If true: print fit summary and chi² const bool drawCanvases = false; // If false: canvases won't be drawn at all -const bool drawCanvases = false; // If false: canvases won't be drawn at all void GainMatchSX3::Begin(TTree * /*tree*/) { @@ -100,13 +99,6 @@ void GainMatchSX3::Begin(TTree * /*tree*/) // frontGain[id][bk][u][d] = gain; // frontGainValid[id][bk][u][d] = true; // } - // 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,6 +127,12 @@ Bool_t GainMatchSX3::Process(Long64_t entry) 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)); @@ -204,52 +202,37 @@ Bool_t GainMatchSX3::Process(Long64_t entry) } } } - for (int i = 0; i < sx3.multi; i++) + + // Only if we found all three channels do we proceed + if (sx3ChUp >= 0 && sx3ChDn >= 0 && sx3ChBk >= 0) { - 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 + // Fill once per correlated set 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] == 4) + // Pick detector ID from one of the correlated hits (all same detector) + int detID = sx3ID[0].first; + + TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", + detID, sx3ChUp, sx3ChDn, sx3ChBk); + TH2F *hist2d = (TH2F *)gDirectory->Get(histName); + if (!hist2d) { - 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 - 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); - } - - hist2d->Fill(sx3EUp + sx3EDn, 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); - } + hist2d = new TH2F(histName, histName, + 400, 0, 16000, 400, 0, 16000); } + + if (sx3EBk > 100 || sx3EUp > 100 || sx3EDn > 100) + { + hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk); + + // Use the correlated triplet directly + dataPoints[{detID, sx3ChBk, sx3ChUp, sx3ChDn}] + .emplace_back(sx3EBk, sx3EUp, sx3EDn); + } + + hist2d->Fill(sx3EUp + sx3EDn, sx3EBk); } } } @@ -260,67 +243,52 @@ Bool_t GainMatchSX3::Process(Long64_t entry) const double GAIN_ACCEPTANCE_THRESHOLD = 0.3; 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}}}}; - std::map upCorrFactor; + double backSlope[MAX_DET][MAX_BK] = {{0}}; + bool backSlopeValid[MAX_DET][MAX_BK] = {{false}}; - // === Gain matching === - std::ofstream outFile("sx3_GainMatchback.txt"); + std::ofstream outFile("sx3_BackGains.txt"); if (!outFile.is_open()) { std::cerr << "Error opening sx3_BackGains.txt for writing!" << std::endl; return; } - // Gain fit using up+dn vs bk - for (const auto &kv : dataPoints) + // === Gain fit: (Up+Dn) vs Back, grouped by [id][bk] === + for (int id = 0; id < MAX_DET; id++) { - // 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; - - std::vector bkE, udE; - - for (const auto &pr : pts) + for (int bk = 0; bk < MAX_BK; bk++) { - double eUp, eDn, eBk; - std::tie(eBk, eUp, eDn) = pr; + std::vector bkE, udE; - if ((eBk < 100) || (eUp < 100) || (eDn < 100)) - continue; // Skip if any energy is less than 100 + // Collect all (Up+Dn, Back) for this id,bk + for (const auto &kv : dataPoints) + { + auto [cid, cbk, u, d] = kv.first; + if (cid != id || cbk != bk) + continue; - bkE.push_back(eBk); - udE.push_back(eUp + eDn); - } + for (const auto &pr : kv.second) + { + double eBk, eUp, eDn; + std::tie(eBk, eUp, eDn) = pr; + if ((eBk < 100) || (eUp < 100) || (eDn < 100)) + continue; - // Fill the TGraph with bkE and udE - // TGraph g(bkE.size(), bkE.data(), udE.data()); - // Fit the graph to a linear function - if (bkE.size() < 5) - continue; // Ensure we have enough points for fitting + bkE.push_back(eBk); + udE.push_back(eUp + eDn); + } + } - const double fixedError = 10.0; // in ADC channels + if (bkE.size() < 5) + continue; // not enough statistics - std::vector xVals, yVals, exVals, eyVals; + // Build graph with errors + const double fixedError = 10.0; // ADC channels + std::vector exVals(udE.size(), 0.0); // no x error + std::vector eyVals(udE.size(), fixedError); // constant y error - // Build data with fixed error - for (size_t i = 0; i < udE.size(); ++i) - { - 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()); + TGraphErrors g(udE.size(), udE.data(), bkE.data(), + exVals.data(), eyVals.data()); TF1 f("f", "[0]*x", 0, 16000); f.SetParameter(0, 1.0); // initial slope @@ -361,39 +329,21 @@ void GainMatchSX3::Terminate() g.Fit(&f, "QNR"); } - 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]) - // { - // if (TMath::Abs(gainArray[id][u][d][bk] - 1) < 0.3) - { - printf("Gain match Det%d Up%dDn%d Backs%d → %.4f \n", id, u, d, bk, gainArray[id][u][d][bk]); - outFile << id << " " << bk << " " << u << " " << d << " " << gainArray[id][u][d][bk] << std::endl; + double slope = f.GetParameter(0); + if (std::abs(slope - 1.0) < 0.3) // sanity check + { + backSlope[id][bk] = slope; + backSlopeValid[id][bk] = true; + outFile << id << " " << bk << " " << slope << "\n"; + printf("Back slope Det%d Bk%d → %.4f\n", id, bk, slope); + } + else + { + std::cerr << "Warning: Bad slope for Det" << id << " Bk" << bk + << " slope=" << slope << 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; - // outFile << id << " " << bk << " " << u << " " << d << " " << gainArray[id][u][d][bk] << std::endl; - // } } - // } - // } - // } - // } - // } outFile.close(); std::cout << "Back gain matching complete." << std::endl; diff --git a/GainMatchSX3Front.C b/GainMatchSX3Front.C index 5513a89..a0417ef 100644 --- a/GainMatchSX3Front.C +++ b/GainMatchSX3Front.C @@ -254,10 +254,10 @@ Bool_t GainMatchSX3Front::Process(Long64_t entry) if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk) && cut1 && cut1->IsInside(sx3EUp / sx3EBk, sx3EDn / sx3EBk)) { - if (backGainValid[sx3.id[i]][sx3ChBk]) - { - sx3EBk *= backGain[sx3.id[i]][sx3ChBk]; - } + // if (backGainValid[sx3.id[i]][sx3ChBk]) + // { + // sx3EBk *= backGain[sx3.id[i]][sx3ChBk]; + // } // Accumulate data for gain matching dataPoints[{sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(sx3EBk, sx3EUp, sx3EDn); }