diff --git a/GainMatch.C b/GainMatch.C index ecafe49..f6ccbb8 100644 --- a/GainMatch.C +++ b/GainMatch.C @@ -7,7 +7,7 @@ #include #include #include - +#include #include #include @@ -26,28 +26,30 @@ PW pw_contr; TVector3 hitPos; bool HitNonZero; +const int MAX_DET = 4; +const int MAX_WEDGE = 16; +const int MAX_RING = 16; + +bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}}; +double gainArray[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}}; + void GainMatch::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); + sx3_contr.ConstructGeo(); pw_contr.ConstructGeo(); + } Bool_t GainMatch::Process(Long64_t entry) { - - // if ( entry > 100 ) return kTRUE; - hitPos.Clear(); HitNonZero = false; - // if( entry > 1) return kTRUE; - // printf("################### ev : %llu \n", entry); - b_sx3Multi->GetEntry(entry); b_sx3ID->GetEntry(entry); b_sx3Ch->GetEntry(entry); @@ -68,41 +70,26 @@ Bool_t GainMatch::Process(Long64_t entry) qqq.CalIndex(); pc.CalIndex(); - // sx3.Print(); - - // ########################################################### Raw data - // //======================= SX3 - // Initialize vector to store pairs of ID and index - std::vector> ID; // first = id, second = index - - // Loop through each entry in sx3.multi + std::vector> ID; for (int i = 0; i < sx3.multi; i++) { - // Store ID and index pair in ID vector ID.push_back(std::pair(sx3.id[i], i)); } - // Check if the ID vector is not empty if (ID.size() > 0) { - // Sort the ID vector by the first element (ID) std::sort(ID.begin(), ID.end(), [](const std::pair &a, const std::pair &b) { return a.first < b.first; }); - // Create a new vector to store IDs that are the same std::vector> sx3ID; - sx3ID.push_back(ID[0]); // Start with the first ID + sx3ID.push_back(ID[0]); bool found = false; - // Loop through the sorted IDs and group by same IDs for (size_t i = 1; i < ID.size(); i++) { if (ID[i].first == sx3ID.back().first) { - // If the ID matches the last one in sx3ID, add to the group sx3ID.push_back(ID[i]); - - // If 3 or more IDs are grouped, set found to true if (sx3ID.size() >= 3) { found = true; @@ -110,7 +97,6 @@ Bool_t GainMatch::Process(Long64_t entry) } else { - // If a new ID is encountered and no group is found, reset the group if (!found) { sx3ID.clear(); @@ -119,21 +105,16 @@ Bool_t GainMatch::Process(Long64_t entry) } } - // If a group of 3 or more IDs is found, process the channels and energies if (found) { int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1; float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0; - // Loop through the grouped IDs for (size_t i = 0; i < sx3ID.size(); i++) { - int index = sx3ID[i].second; // Get the index from the pair - - // Categorize channel and energy + int index = sx3ID[i].second; if (sx3.ch[index] < 8) { - // If channel is less than 8, assign it to up or down if (sx3.ch[index] % 2 == 0) { sx3ChDn = sx3.ch[index]; @@ -147,88 +128,99 @@ Bool_t GainMatch::Process(Long64_t entry) } else { - // If channel is >= 8, assign it to back sx3ChBk = sx3.ch[index]; sx3EBk = sx3.e[index]; } } - // make a histogram for Sx3 back energy vs sum of sx3 front energies hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk); - // Further energy matching calculations can be added here... } } - // //======================= QQQ - TH2F *hist1 = NULL; - // Loop through each entry in qqq.multi + for (int i = 0; i < qqq.multi; i++) { 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; float eWedge = 0.0; float eRing = 0.0; - if (qqq.ch[i] < qqq.ch[j]) + if (qqq.ch[i] < 16 && qqq.ch[j] >= 16) { - chRing = qqq.ch[j] - 16; - eRing = qqq.e[j]; chWedge = qqq.ch[i]; eWedge = qqq.e[i]; + chRing = qqq.ch[j] - 16; + eRing = qqq.e[j]; + } + else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16) + { + chWedge = qqq.ch[j]; + eWedge = qqq.e[j]; + chRing = qqq.ch[i] - 16; + eRing = qqq.e[i]; } else - { - chRing = qqq.ch[i]; - eRing = qqq.e[i]; - chWedge = qqq.ch[j] - 16; - eWedge = qqq.e[j]; - } - // printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing); - hQQQFVB->Fill(qqq.e[i], qqq.e[j]); - // 1. Create/get individual 2D histogram - TString histName = Form("hQQQFVB_r%d_w%d_id%d", chRing, chWedge, qqq.id[i]); + continue; + + hQQQFVB->Fill(eWedge, eRing); + + TString histName = Form("hQQQFVB_id%d_r%d_w%d", qqq.id[i], chRing, chWedge); TH2F *hist2d = (TH2F *)gDirectory->Get(histName); if (!hist2d) { - hist2d = new TH2F(histName, Form("QQQ GainMatch ID%d R%d W%d;Wedge E;Ring E", qqq.id[i], chRing, chWedge), - 400, 0, 16000, 400, 0, 16000); + hist2d = new TH2F(histName, Form("QQQ GainMatch Det%d R%d W%d;Wedge E;Ring E", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000); } hist2d->Fill(eWedge, eRing); } } } + return kTRUE; } + void GainMatch::Terminate() { + + TFile *cutFile = TFile::Open("qqqcorr.root"); + TCutG *cut = (TCutG *)cutFile->Get("qqqcorr"); + const int MAX_DET = 4; + const int MAX_RING = 16; + const int MAX_WEDGE = 16; + + // Store gains and validity + static double gainArray[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}}; + static bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}}; + + std::ofstream outFile("qqq_gainmatch.txt"); + if (!outFile.is_open()) + { + printf("Error opening output file!"); + return; + } + + // Collect all (wedgeE, ringE) pairs per detector/ring/wedge + std::map, std::vector>> dataPoints; TIter next(gDirectory->GetList()); TObject *obj; - while ((obj = next())) { if (!obj->InheritsFrom("TH2F")) continue; - - TH2F *hist2d = (TH2F *)obj; + TH2F *hist2d = static_cast(obj); TString name = hist2d->GetName(); - if (!name.BeginsWith("hQQQFVB_r")) + if (!name.BeginsWith("hQQQFVB_id")) + continue; + if (hist2d->GetEntries() < 100) continue; - if (hist2d->GetEntries() < 1000) - continue; - - std::vector wedge_vals; - std::vector peak_vals; + int id, ring, wedge; + sscanf(name.Data(), "hQQQFVB_id%d_r%d_w%d", &id, &ring, &wedge); for (int binX = 1; binX <= hist2d->GetNbinsX(); ++binX) { - TH1D *projY = hist2d->ProjectionY("_py", binX, binX); + TH1D *projY = hist2d->ProjectionY("_py", binX, binX + 1); if (projY->GetEntries() < 30) { delete projY; @@ -240,34 +232,96 @@ void GainMatch::Terminate() projY->Fit(fit, "QNR"); double wedgeE = hist2d->GetXaxis()->GetBinCenter(binX); - double ringPeak = fit->GetParameter(1); + double ringE = fit->GetParameter(1); + dataPoints[{id, ring, wedge}].emplace_back(wedgeE, ringE); - wedge_vals.push_back(wedgeE); - peak_vals.push_back(ringPeak); - - delete projY; delete fit; + delete projY; + } + } + + // Fit slopes with sigma-clipping + for (auto &kv : dataPoints) + { + auto [id, ring, wedge] = kv.first; + auto &pts = kv.second; + if (pts.size() < 5) + continue; + + // Build vectors + std::vector wE, rE; + for (auto &pr : pts) + { + wE.push_back(pr.first); + rE.push_back(pr.second); } - if (wedge_vals.size() >= 5) + // Initial fit + TGraph *g0 = new TGraph(wE.size(), wE.data(), rE.data()); + TF1 *f0 = new TF1("f0", "[0]*x + [1]", 0, 16000); + g0->Fit(f0, "QNR"); + double m0 = f0->GetParameter(0), b0 = f0->GetParameter(1); + + // Clip to cut + std::vector wC, rC; + for (size_t i = 0; i < wE.size(); ++i) { - TGraph *gr = new TGraph(wedge_vals.size(), &wedge_vals[0], &peak_vals[0]); - gr->SetName(name + "_fit"); - TF1 *line = new TF1("line", "pol1", 0, 16000); - gr->Fit(line, "Q"); + if (cut->IsInside(wE[i], rE[i])) + { + wC.push_back(wE[i]); + rC.push_back(rE[i]); + } + } - double gain = line->GetParameter(1); - double offset = line->GetParameter(0); - printf("Gain match %s → Gain = %.4f, Offset = %.2f\n", name.Data(), gain, offset); + // Final fit + if (wC.size() >= 5) + { + TGraph *g1 = new TGraph(wC.size(), wC.data(), rC.data()); + TF1 *f1 = new TF1("f1", "[0]*x + [1]", 0, 16000); + g1->Fit(f1, "QNR"); + gainArray[id][ring][wedge] = f1->GetParameter(0); + gainValid[id][ring][wedge] = true; + delete g1; + delete f1; + } + delete g0; + delete f0; + } - TCanvas *c1 = new TCanvas(name + "_c", name + "_c"); - gr->SetTitle(Form("Gain Match: %s", name.Data())); - gr->GetXaxis()->SetTitle("Wedge Energy"); - gr->GetYaxis()->SetTitle("Ring Energy"); - gr->Draw("AP"); - line->Draw("same"); - c1->SaveAs(Form("%s_fit.png", name.Data())); - delete c1; + // Write out valid entries + for (int id = 0; id < MAX_DET; ++id) + { + for (int ring = 0; ring < MAX_RING; ++ring) + { + for (int wedge = 0; wedge < MAX_WEDGE; ++wedge) + { + if (!gainValid[id][ring][wedge]) + continue; + outFile << id << " " << wedge << " " << ring << " " << gainArray[id][ring][wedge] << std::endl; + // printf("Gain match Det%d Ring%d Wedge%d → %.4f \n", id, ring, wedge, gainArray[id][ring][wedge]); + } + } + } + outFile.close(); + printf("Gain matching complete (with sigma-clipping)"); + + // === Plot all gain-matched QQQ points together with a 2D histogram === + TH2F *hAll = new TH2F("hAll", "All QQQ Gain-Matched;Corrected Wedge E;Ring E", + 400, 0, 16000, 400, 0, 16000); + + // Fill the combined TH2F with corrected data + for (auto &kv : dataPoints) + { + int id, ring, wedge; + std::tie(id, ring, wedge) = kv.first; + if (!gainValid[id][ring][wedge]) + continue; + auto &pts = kv.second; + for (auto &pr : pts) + { + double corrWedge = pr.first * gainArray[id][ring][wedge]; + double ringE = pr.second; + hAll->Fill(corrWedge, ringE); } } }