diff --git a/Calibration.C b/Calibration.C index 2315def..4e28179 100644 --- a/Calibration.C +++ b/Calibration.C @@ -97,9 +97,9 @@ void Calibration::Begin(TTree * /*tree*/) } else { - int id, bk, u, d; + int id, bk; double gain; - while (infile >> id >> bk >> u >> d >> gain) + while (infile >> id >> bk >> gain) { backGain[id][bk] = gain; backGainValid[id][bk] = (gain > 0); @@ -253,71 +253,69 @@ Bool_t Calibration::Process(Long64_t entry) int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1; float sx3EUp = 0.0f, sx3EDn = 0.0f, sx3EBk = 0.0f; - for (auto &p : sx3ID) + for (size_t i = 0; i < sx3ID.size(); i++) { - int index = p.second; - int ch = sx3.ch[index]; - float e = sx3.e[index]; + int index = sx3ID[i].second; - if (ch < 8) + if (sx3.ch[index] < 8) { - if ((ch % 2) == 0) // even -> down + if ((sx3.ch[index] % 2) == 0) // even -> down { - sx3ChDn = ch; - sx3EDn = e; + sx3ChDn = sx3.ch[index]; + sx3EDn = sx3.e[index]; } else // odd -> up { - sx3ChUp = ch; - sx3EUp = e; + sx3ChUp = sx3.ch[index]; + sx3EUp = sx3.e[index]; } } else { - sx3ChBk = ch; - sx3EBk = e; - } - } - - bool haveFrontPair = (sx3ChUp >= 0 || sx3ChDn >= 0); - bool haveBack = (sx3ChBk >= 0); - double GM_EUp = 0.0, GM_EDn = 0.0, calibEBack = 0.0; - - if (haveBack) - { - // --- 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[0].first][sx3ChBk][sx3ChUp][sx3ChDn]->Fill(sx3EUp); - - // --- If gain is available, also fill calibrated energy - if (frontGainValid[sx3ID[0].first][sx3ChBk][sx3ChUp][sx3ChDn]) - { - GM_EUp = frontGain[sx3ID[0].first][sx3ChBk][sx3ChUp][sx3ChDn] * sx3EUp; - if (GM_EUp > 50.0) - hSX3Spectra[sx3ID[0].first][sx3ChBk][sx3ChUp][sx3ChDn]->Fill(GM_EUp); // optional: mixes raw+calib - } - // --- If back gain is available, also fill calibrated energy - hsx3E_raw->Fill(sx3EBk); - - if (backGainValid[sx3ID[0].first][sx3ChBk]) - { - calibEBack = backGain[sx3ID[0].first][sx3ChBk] * sx3EBk; - if (calibEBack > 50.0) - hsx3E_calib->Fill(calibEBack); // optional: mixes raw+calib + sx3ChBk = sx3.ch[index]; + sx3EBk = sx3.e[index]; } - // 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); + bool haveFrontPair = (sx3ChUp >= 0 || sx3ChDn >= 0); + bool haveBack = (sx3ChBk >= 0); + double GM_EUp = 0.0, GM_EDn = 0.0, calibEBack = 0.0; - if (GM_EUp > 50.0 && sx3EBk > 50.0) + if (haveBack) { - sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, GM_EUp, sx3EDn); - hitPos = sx3_contr.GetHitPos(); - HitNonZero = true; + // --- 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[i].first][sx3ChBk][sx3ChUp][sx3ChDn]->Fill(sx3EUp); + + // --- If gain is available, also fill calibrated energy + if (frontGainValid[sx3ID[i].first][sx3ChBk][sx3ChUp][sx3ChDn]) + { + GM_EUp = frontGain[sx3ID[i].first][sx3ChBk][sx3ChUp][sx3ChDn] * sx3EUp; + if (GM_EUp > 50.0) + hSX3Spectra[sx3ID[i].first][sx3ChBk][sx3ChUp][sx3ChDn]->Fill(GM_EUp); // optional: mixes raw+calib + } + // --- If back gain is available, also fill calibrated energy + hsx3E_raw->Fill(sx3EBk); + + if (backGainValid[sx3ID[i].first][sx3ChBk]) + { + calibEBack = backGain[sx3ID[i].first][sx3ChBk] * sx3EBk; + if (calibEBack > 50.0) + hsx3E_calib->Fill(calibEBack); // optional: mixes raw+calib + } + + // Keep the other diagnostic plots + hsx3IndexVE_gm->Fill(sx3.index[sx3ID[i].second], GM_EUp); + hSX3->Fill(sx3ChDn + 4, sx3ChBk); + hSX3->Fill(sx3ChUp, sx3ChBk); + hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk); + + if (GM_EUp > 50.0 && sx3EBk > 50.0) + { + sx3_contr.CalSX3Pos(sx3ID[i].first, sx3ChUp, sx3ChDn, sx3ChBk, GM_EUp, sx3EDn); + hitPos = sx3_contr.GetHitPos(); + HitNonZero = true; + } } } } diff --git a/GainMatchSX3.C b/GainMatchSX3.C index 483cbf3..3f58e2c 100644 --- a/GainMatchSX3.C +++ b/GainMatchSX3.C @@ -283,7 +283,7 @@ void GainMatchSX3::Terminate() continue; // not enough statistics // Build graph with errors - const double fixedError = 10.0; // ADC channels + const double fixedError = 0.0; // ADC channels std::vector exVals(udE.size(), 0.0); // no x error std::vector eyVals(udE.size(), fixedError); // constant y error @@ -291,7 +291,7 @@ void GainMatchSX3::Terminate() exVals.data(), eyVals.data()); TF1 f("f", "[0]*x", 0, 16000); - f.SetParameter(0, 1.0); // initial slope + // f.SetParameter(0, 1.0); // initial slope if (drawCanvases) { @@ -329,7 +329,7 @@ void GainMatchSX3::Terminate() g.Fit(&f, "QNR"); } - double slope = f.GetParameter(0); + double slope = 1/f.GetParameter(0); if (std::abs(slope - 1.0) < 0.3) // sanity check { backSlope[id][bk] = slope; diff --git a/GainMatchSX3Front.C b/GainMatchSX3Front.C index a0417ef..61de72e 100644 --- a/GainMatchSX3Front.C +++ b/GainMatchSX3Front.C @@ -44,17 +44,15 @@ bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; double uvdslope[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; // ==== Configuration Flags ==== const bool interactiveMode = false; // If true: show canvas + wait for user -const bool interactiveMode1 = 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 drawCanvases1 = false; // If false: canvases won't be drawn at all 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); + hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 800, 0, 16000, 800, 0, 16000); + hSX3FvsB_g = new TH2F("hSX3FvsB_g", "SX3 Front vs Back; Front E; Back E", 800, 0, 16000, 800, 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); @@ -84,27 +82,6 @@ void GainMatchSX3Front::Begin(TTree * /*tree*/) 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; - // if (backGain[id][bk][u][d] > 0) - // backGainValid[id][bk][u][d] = true; - // else - // backGainValid[id][bk][u][d] = false; - // } - // infile.close(); - // std::cout << "Loaded back gains from " << filename << std::endl; std::string filename = "sx3_BackGains.txt"; @@ -115,7 +92,7 @@ void GainMatchSX3Front::Begin(TTree * /*tree*/) return; } - int id, bk, u, d; + int id, bk; double gain; while (infile >> id >> bk >> gain) { @@ -254,10 +231,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); } @@ -290,7 +267,7 @@ void GainMatchSX3Front::Terminate() auto [id, bk, u, d] = kv.first; const auto &pts = kv.second; - if (pts.size() < 500) + if (pts.size() < 50) continue; std::vector uE, dE, udE, corrBkE; @@ -322,20 +299,20 @@ void GainMatchSX3Front::Terminate() // 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 + double x = uE[i]; // front energy + double y = dE[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 + exVals.push_back(fixedError); // error in up energy + // eyVals.push_back(fixedError); // error in down 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 + f.SetParameter(0, -1.0); // Initial guess if (drawCanvases) { @@ -353,9 +330,9 @@ void GainMatchSX3Front::Terminate() 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; + 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; } if (interactiveMode) @@ -380,87 +357,6 @@ void GainMatchSX3Front::Terminate() // printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]); } - for (const auto &kv : dataPoints) - { - auto [id, bk, u, d] = kv.first; - const auto &pts = kv.second; - std::vector uvals, dvals; - double front = frontGain[id][bk][u][d]; - - if (pts.size() < 5) - continue; - - 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 - eUp *= front; - // eDn *= frontGain[id][bk][u][d]; - - uE.push_back(eUp / eBkCorr); - dE.push_back(eDn / eBkCorr); - corrBkE.push_back(eBkCorr); - } - for (size_t i = 0; i < uE.size(); ++i) - { - uvals.push_back(uE[i]); - dvals.push_back(dE[i]); - } - - TGraph g1(uvals.size(), uvals.data(), dvals.data()); - g1.SetMarkerStyle(20); - g1.SetMarkerColor(kBlue); - g1.SetTitle(Form("Det %d: U%d D%d B%d", id, u, d, bk)); - - // Fit function (straight line through origin) - TF1 f1("f1", "[0]*x+[1]", 0, 16000); - f1.SetParameter(0, -1.0); - - if (drawCanvases1) - { - TCanvas *c1 = new TCanvas(Form("c_%d_%d_%d_%d", id, bk, u, d), "Fit", 800, 600); - g1.SetTitle(Form("Detector %d: U%d D%d B%d", id, u, d, bk)); - g1.SetMarkerStyle(20); - g1.SetMarkerColor(kBlue); - g1.Draw("AP"); - - g1.Fit(&f1, interactiveMode1 ? "Q" : "QNR"); // 'R' avoids refit, 'N' skips drawing - - if (verboseFit) - { - double chi2 = f1.GetChisquare(); - int ndf1 = f1.GetNDF(); - double reducedChi2 = (ndf1 != 0) ? chi2 / ndf1 : -1; - - std::cout << Form("Det%d U%d D%d B%d → Gain: %.4f | χ²/ndf = %.2f/%d = %.2f", - id, u, d, bk, f1.GetParameter(0), chi2, ndf1, reducedChi2) - << std::endl; - } - - if (interactiveMode1) - { - c1->Update(); - gPad->WaitPrimitive(); - } - else - { - c1->Close(); // Optionally avoid clutter in batch - } - } - else - { - g1.Fit(&f1, "QNR"); - } - - // Save slope - uvdslope[id][bk][u][d] = f1.GetParameter(0); - printf("UvD slope Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, uvdslope[id][bk][u][d]); - } - outFile.close(); std::cout << "Gain matching complete." << std::endl;